A  FINITE  ELEMENT  MODEL  OF 
FLUID  FLOW  IN  JOINTED  ROCK 


BY 

MARK  A.  JAMES 
B.S.,  Kansas  State  University,  1986 


A  THESIS 

submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree 

MASTER  OF  SCIENCE 


Department  of  Mechanical  Engineering 

KANSAS  STATE  UNIVERSITY 
Manhattan,  Kansas 

1988 


Approved  by: 
Major  Professor 


A1120S  2DS1EM 

Acknowledgments 

.  / 

Daniel  Swenson,  my  advisor,  deserves  more  than  my  thanks  for  his 

patience  and  for  his  moral  and  technical  support.   He  taught  me  much 

more  than  finite  elements  and  computer  programming. 

Shun  Lung  Su  developed  the  fluid  model  which  was  implemented  in 
DRACULA. 

Thanks  should  go  to  Don  Brown  for  working  to  get  the  financial 
support  to  finish  this  work  and  for  helping  me  understand,  at  least  to  some 
degree,  what  I  am  modelling. 

My  parents  always  support  me  in  everything  I  do.   Thank  you. 


M 


Table  of  Contents 

Table  of  Contents  iii 

List  of  Figures  v 

List  of  Tables  viii 

Chapter  I         1 

1.1  The  HDR  Concept  1 

1.2  The  HDR  Program   2 

1.3  Previous  work  4 

1.4  Objectives  and  Scope  6 

Chapter  II       11 

2.1  Element  Derivations  12 

2.1.1  Elasto-Dynamic  Structural  Element  12 

2.1.2  Fluid  Flow  Element  15 

2.1 .3  Interface  and  Boundary  Spring  Elements   18 

2.2  Joint  Model  21 

2.3  Solution  Method  22 

2.4  Coupling  between  Fluid  and  Structural  Models 26 

2.5  Tracer  Model  27 

Chapter  III     33 

3.1  Implementation    33 

3.2  Modelling  Approach  37 

3.2.1  Joint  Model  Implementation 38 

3.2.2  Reservoir  Boundary  Conditions  39 

Chapter  IV      44 


ii 


4.1  Steady  State  Verification  Problem  44 

4.2  Transient  Verification  Problem 46 

4.3  Summary  49 

Chapter  V  67 

5.1  Low  Pressure  Steady  State  Problem  67 

5.2  High  Pressure  Steady  State  Problem  69 

5.3  Transient  Flow  Problem  72 

5.4  Summary  74 

Chapter  VI  93 

6.1  Summary  93 

6.2  Conclusions  95 

6.3  Recommendations  97 

References  100 


IV 


List  of  Figures 

1.1  Hot  Dry  Rock  Concept  8 

1.2  Anticipated  Phase  II  Fracture  System  9 

1.3  Cubic  Flow  Law  10 

2.1  Differential  Element  for  Conservation  of  Mass   28 

2.2  Typical  Joint  Opening  Law  Traction  -  Displacement 
Relation  29 

2.3  Typical  Boundary  Spring  Traction  —  Displacement 

Relation  29 

2.4  "Bed  of  Nails"  Joint  Model 30 

2.5  Hydraulic  Aperture  and  Tracer  Aperture 31 

3.1  DRACULA  Menu  Concept  40 

3.2  Joint  Opening  Law,  Interface  Materials  One  and  Two  41 

3.3  Application  Problems  Finite  Element  Mesh  42 

3.4  Boundary  Springs,  Interface  Materials  Three  and  Four....  43 

4.1  Verification  Problems  Finite  Element  Mesh 52 

4.2  Verification  Problems  Boundary  Springs 53 

4.3a     Far  Field  Flow  Loss  Boundary  Condition 54 

4.3b     Flow  Rate  Boundary  Condition 54 

4.4  Flow  Rate  Contour  Plot  55 

4.5  Pressure  Contour  Plot  56 

4.6  Joint  Effective  Opening  Stress  56 

4.7  Y  Stress  Line  Plot  57 


4.8  Displaced  Mesh  Plot  58 

4.9  Y  Stress  Line  Plot,  time  =  0.0  (transient  verification 
problem)  59 

4.10  Time  History  Plot,  Initial  9  MPa  Solution 60 

4.11  Time  History  Plot  of  Displaced  Mesh  Plots,  t-step  =  0.125...  61 

4.12  Initial  Displaced  Mesh  Plot  63 

4.13  This  figure  intentionally  left  blank  64 

4.14  Flow  Rate  and  Pressure  Transient  Solution  65 

4.15  Fluid  Flow  Rate  Solution  for  Different  values  of  t-step 66 

5.1  Boundary  Conditions  (low  pressure  steady  state 

problem)  76 

5.2  Pressure  Contour  Plot  77 

5.3  Flow  Rate  Contour  Plot  78 

5.4  Flow  Rate  Contour,  Zoom  Around  Extraction  Well 78 

5.5  Displaced  Mesh  Plot 79 

5.6  Displacement  History  Plots  80 

5.7  Line  Plot,  Y  Stress  83 

5.8  Pressure  Contour  Plot  (high  pressure  steady  state 

problem)  84 

5.9  Y  Stress  Line  Plot  85 

5.10  X  Stress  Line  Plot 86 

5.11  Displacement  History  Plots  87 

5.12  Experimental  Extraction  Well  Pressure  Results, 
Experiment  #2070  88 

5.13  Applied  Pressure  Boundary  Condition  (transient 
application  problem) 89 


VI 


5.14  Pressure  Contour  Plot,  First  Transient  Step 90 

5.15  Comparison  of  Calculated  and  Experimental  Output 
Pressure  Histories  91 


VII 


list  of  Tables 

4.1  Reservoir  Properties 50 

4.2  Joint  Volume  Change 51 


VIM 


Chapter  I 
Introduction 

In  the  early  1970's  the  United  States  experienced  its  first  energy- 
crisis.   At  that  time  about  97%  of  the  world's  energy  needs  were  being 
supplied  by  fossil  fuels.  The  need  for  the  United  States  to  become  less 
dependent  on  petroleum  sources  led  to  an  extensive  search  for  alternate 
energy  sources.   One  of  these  alternatives  is  geothermal  energy. 

Geothermal  energy  is  the  natural  heat  produced  within  the  Earth's 
crust  by  slowly  decaying,  naturally  occurring,  radioactive  isotopes  (i.e. 
uranium,  thorium  and  potassium)  and  by  conduction  from  the  hotter 
interior  regions.    This  energy  is  stored  in  several  forms:  hydrothermal 
reservoirs  (naturally-occurring  pockets  of  steam  or  hot  water), 
geopressured  reservoirs,  magma  reservoirs  and  hot  dry  rock  reservoirs. 
This  thesis  presents  a  computer  model  suitable  for  analyzing  hot  dry  rock 
reservoirs  in  two  dimensions  using  the  finite  element  method. 

1.1       The  HDR  Concept 

In  a  geothermal  energy  context,  hot  dry  rock  (HDR)  is  defined  as 
naturally  heated,  unmelted  crustal  rock,  that  does  not  produce  natural 
steam  or  hot  water  at  commercially  useful  rates.   HDR  exists  everywhere 
beneath  the  earth's  surface;  however,  the  quality  of  the  resource  depends  on 
the  local  temperature  gradient.  At  a  given  location,   the  temperatures  at 
economical  depths  may  not  be  high  enough  for  electric  power  generation, 
but  would  almost  everywhere  be  suitable  for  direct  use  in  agriculture  or 
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food  and  chemical  processing,  or  for  supplemental  energy  generation  or 
space  heating  (Armstead,  1978). 

The  energy  is  extracted  from  a  region  of  hot  rock  by  circulating 
pressurized  water  in  a  closed  loop  through  a  man-made  fracture  system 
created  by  hydraulically  fracturing  the  rock  between  two  well-bores.   The 
useful  heat  is  recovered  at  the  surface  through  heat  exchangers;  the  cooled 
water  is  reinjected  to  recirculate  through  the  underground  loop  (see  Figure 
1.1).  Current  estimates  show  that  at  depths  <10  km,  the  HDR  reservoir  base 
contains  -32  million  Quads  of  energy.1   Two  percent  of  this  energy  would 
supply  the  United  States'  nontransportation  needs  for  about  2000  years  at 
present  consumption  rates  (Los  Alamos  annual  report,  1979). 

1.2       The  HDR  Program 

The  HDR  program  originated  at  Los  Alamos  National  Laboratory 
(LANL)  in  1970  as  a  feasibility  study  of  underground  heat  extraction  from 
low  permeability  rock.  The  work  at  LANL  continued,  and  in  1977  the 
world's  first  HDR  geothermal  energy  system  (Phase  I)  was  completed  at 
Fenton  Hill,  New  Mexico.   In  1979  the  system  was  enlarged  using  hydraulic 
fracturing  and  a  series  of  additional  flow  tests  were  run.  The  tests  led  to  a 
reasonably  well  defined  model  for  the  system  (Los  Alamos  annual  report, 
1981).  During  flow  tests  in  1980  and  1981  heat  extraction  for  the  Phase  I 
system  rose  to  ~5  MWt  and  the  rate  of  water  loss  decreased  to  about  1.5%  of 
the  flow  rate.  An  electric  generator  was  operated  at  its  full  rated  capacity  of 


1  1  Quad  =  1  quadrillion  (1015)  Btu  =  334  MW  centuries  ~1018  J. 


60  kwA  and  exceeded  its  overall  design  efficiency  of  5.7%. 

During  1979  a  larger,  deeper,  hotter  Phase  II  system  was  begun.   The 
initial  drilling  for  the  two  Phase  II  wells  was  in  1980.   Initial 
hydrofracturing  from  the  injection  well  in  1983  produced  a  fracture  system 
quite  different  than  observed  in  the  Phase  I  wells.   Figure  1 .2  shows  the 
anticipated  Phase  II  fracture  system  based  on  an  analysis  of  the  Phase  I 
and  other  HDR  systems.   Two  types  of  joints  make  up  the  reservoir  flow 
paths:  vertical  shear  joints  that  are  initially  closed  and  require  high 
pressures  to  open,  and  inclined  tensile  joints  that  are  also  initially  closed, 
but  require  lower  pressures  to  open.   The  actual  Phase  II  fractured  system 
determined  from  the  location  of  microseismic  events  is  three  dimensional 
rather  than  planar,  inclined  rather  than  vertical,  and  did  not  hydraulically 
connect  the  two  wells.    Rather  than  continue  with  more  hydrofracturing, 
the  injection  well  was  directionally  redrilled  in  1985  to  intercept  the  fracture 
system  created  from  the  production  well.    This  redrilling  was  successful, 
and  subsequent  hydrofracturing  improved  hydraulic  communication 
between  the  two  wells  (Franke,  1988). 

In  May  and  June  1986  the  Initial  Closed  loop  Flow  Test  (ICFT)  was 
run  on  the  Phase  II  system  to  obtain  operating  characteristics  needed  to 
plan  the  year  long  Long  Term  Flow  Test  (LTFT).  The  LTFT  will  be  used  to 
predict  thermal  draw-down  as  well  as  other  long  term  effects  required  for 
thorough  reservoir  evaluation.   The  30  day  ICFT  succeeded  with  final 
production  of  about  10  MWt  at  192°  C,  while  injecting  285  gpm  at  4600  psi. 

The  final  water  loss  rate  and  flow  impedance  were  high,  27%  and  18 


psi/gpm  respectively,  but  were  still  declining  (Kelkar  et  al.,  1987).   Birdsell 
and  Robinson  (1988)  have  had  moderate  success  modeling  the  reservoir  as 
an  equivalent  porous  medium.  They  used  the  code  FEHM  developed  by 
Zyvolosky  et  al.  (1988)  which  includes  heat  and  mass  transport  effects  but 
currently  does  not  include  any  deformable  rock  effects.   FEHM  also  has  the 
capability  to  perform  uncoupled  tracer  calculations. 

One  of  the  primary  goals  of  the  HDR  program  has  been  to 
characterize  the  fractured  reservoirs  to  allow  economic  analysis  and 
prediction  of  reservoir  production,  and  ultimately  to  determine  the 
economic  feasibility  of  a  commercial  HDR  energy  production  site. 
Characterization  of  an  HDR  reservoir  usually  proceeds  by  analyzing 
seismic,  temperature,  flow  rate,  and  tracer  data.   Many  new  diagnostic 
methods  have  been  developed  from  the  HDR  programs  in  the  U.S.  and  other 
countries  including  statistical  microseismic  event  analysis,  porous  flow 
models,  tracer  analysis,  as  well  as  steady  state  and  transient  discrete  crack 
flow  models.   This  thesis  presents  a  method  for  analyzing  transient  fluid 
flow  in  discrete  cracks  with  coupling  between  the  fluid  flow  and  rock 
deformations.   In  this  model,  both  the  crack  porosity  and  flow  conductance 
are  nonlinear  functions  of  the  fluid  pressure. 

1.3       Previous  work 

The  traditional  approach  to  modeling  fluid  flow  through  discrete 
fractures  has  been  to  assume  viscous,  incompressible  flow  between  smooth 
parallel  plates  (Snow,  1965).   This  has  come  to  be  known  as  the  "cubic  law," 
in  which  the  volumetric  flow  rate  is  proportional  to  the  pressure  gradient 


and  the  joint  aperture  cubed  (see  Figure  1.3).  Witherspoon  et  al.  (1980)  and 
Ryan  (1987)  verified  the  validity  of  the  "cubic  law"  in  laboratory  work  for 
laminar  flow  between  parallel  planar  plates.   The  cubic  law  was  verified  for 
open  joints  and  for  closed  joints  down  to  a  minimum  of  4  urn.    Witherspoon 
et  al.  (1980)  also  discussed  deviations  from  the  parallel  plate  model  and 
recommends  using  a  factor  of  roughness  in  the  flow  equation,  particularly 
for  high  (>  10  MPa)  crack  closure  stresses.  Values  for  the  factor  of 
roughness  vary  from  1.04  to  1.78.  Brown  (1987)  also  recommends  this.   Su 
(1988)  developed  a  finite  element  model  for  1-D  fluid  flow  based  on  the  cubic 
law  and  included  an  uncoupled  tracer  model. 

Several  finite  element  models  of  coupled  fluid  flow  in  fractured  rock 
masses  have  been  developed  recently.  Noorishad  et  al.  (1982)  developed  a 
code  to  solve  two-dimensional  quasi-static  saturated  porous  media.   Their 
results  show  that  significant  differences  in  pressure  distribution  and  flow 
rate  occur  due  to  coupling  fluid  flow  with  rock  deformations  compared  with 
the  uncoupled  flow  solution.   Hilber  and  Taylor  (1976)  developed  a  dynamic 
code  for  discrete  fractures  that  takes  into  account  the  inertial  effects  of  both 
the  fluid  and  fracture  movements.   The  code  has  been  used  to  study  seismic 
events  along  predetermined  faults  due  to  fluid  injection.   Cundall  (1982) 
developed  the  Fluid  Rock  Interaction  Program  (FRIP)  to  analyze  the 
dynamic  behavior  of  discrete  flow  paths  and  discrete  blocks.   The  blocks 
interact  with  each  other  as  well  as  with  the  fluid.   FRIP  has  been  used  for 
several  transient  response  studies  of  a  geothermal  reservoir  (Pine  and 
Batchelor,  1984),  (Pine  and  Ledingham,  1984),  (Pine  and  Cundall,  1985). 


Asgian  (1988  a)  also  developed  a  discrete  flow  model.  The  FFFLOW 
model  solves  for  quasi-static  rock  deformations  and  transient  fluid  flow  in 
two  dimensions.    The  continuum  rock  masses  are  linear  elastic;  the  joints 
are  nonlinear  elastic  and  include  both  normal  and  shear  rock  stresses  as 
well  as  fluid  pressures.   The  results  indicate  that  the  enhanced 
permeability  zone  is  not  the  same  as  the  pressurized  zone.  Asgian  (1988  b) 
studied  the  transient  response  due  to  different  pumping  rates.   The  studies 
indicate  that  the  peak  response  (peak  pressure,  slippage,  and  aperture 
change)  is  less  intense  for  lower  pumping  rates  than  for  higher  pumping 
rates  of  equal  volumes  of  fluid. 

1.4       Objectives  and  Scope 

The  primary  objective  of  this  thesis  is  to  develop  a  finite  element 
model  of  fluid  flow  through  fractured  rock.   The  finite  element  method  was 
chosen  because  the  governing  equations  for  both  the  rock  masses  and  fluid 
flow  are  well  understood,  and  implementation  is  straight  forward.   The 
model  developed  is  capable  of  solving  the  highly  nonlinear  equations  and  is 
capable  of  solving  very  large  problems.  Extensive  use  of  the  interactive  pre- 
and  post-processing  developed  by  Swenson  (1985)  makes  this  model  a  finite 
element  "analysis  system."   This  thesis  extends  Su's  (1988)  work  to  include 
the  coupled  fluid  flow  -  rock  deformation  model.   This  model  is  the  second 
step  in  developing  a  coupled  fluid  flow  -  rock  deformation  -  heat  transfer 
model  capable  of  modeling  HDR  reservoirs  in  two  dimensions. 

The  HDR  reservoir  is  modeled  as  a  horizontal  plane  with  discrete 
flow  paths  to  model  the  fluid  flow.   The  jointed  rock  mass  forming  the  flow 


paths  are  not  porous  media.   Initially  open  joints  can  be  modeled  as  empty 
fluid  elements  that  fill  up  with  time.   Because  of  the  robust  solution 
algorithm,  free  floating  rock  masses  can  effectively  be  modeled. 

In  this  model,  as  in  the  model  developed  by  Su  (1988)  the  fluid  density 
and  viscosity  are  assumed  constant.    In  the  next  model,  which  will  include 
heat  transfer  in  both  the  structure  as  well  as  the  fluid,  the  temperature 
dependent  properties  of  the  fluid  will  be  handled  more  generally. 

The  remainder  of  the  thesis  is  divided  into  6  chapters.   Chapter  II 
presents  the  finite  element  derivations  for  the  rock  and  fluid  models,  and 
Chapter  III  explains  how  these  two  independent  models  are  coupled. 
Chapter  IV  describes  the  problems  used  to  verify  the  model,  and  Chapter  V 
explains  results  of  problems  that  model  the  HDR  reservoir  at  Fenton  Hill, 
New  Mexico.    Chapter  VI  gives  a  summary  of  the  thesis  and  conclusions 
about  the  results  presented. 
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Figure  1 .1 :       Hot  Dry  Rock  Concept 

(From  LANL  FY83  Report) 
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Figure  1.2:       Anticipated  Phase  II  Fracture  Syst 
(From  LANL  FY83  Report) 
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Figure  1.3: 


Cubic  Flow  Law 
(From  Asgian,  1988a) 
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Chapter  II 
Finite  Element  Model 

A  finite  element  model  of  a  problem  gives  a  piecewise  approximation 
to  the  governing  equations.   The  region  in  which  the  solution  is  desired  is 
divided  into  discrete  elements  and  an  approximate  solution  is  assumed  over 
the  discrete  element  region.   The  contribution  of  each  element  is  then  added 
to  a  global  system  matrix,  which  can  be  solved  for  the  nodal  unknowns. 
The  first  section  of  this  chapter  develops  the  finite  element  method  (FEM) 
for  each  of  the  four  element  types  used  in  this  thesis. 

A  joint  model  must  have  certain  characteristics  to  model  the  complex 
relationship  between  the  fluid  flow  and  the  joint  opening.   The  second  part 
of  this  chapter  discusses  some  rock  mechanics  fundamentals  and  presents 
the  joint  model  implementation. 

In  this  application  of  FEM  we  are  interested  in  solving  for  a  set  of 
nodal  unknowns  (displacements  and  pressures)  that  interact  with  each 
other  through  a  highly  nonlinear  relationship.    The  joint  displacements  are 
a  function  of  joint  fluid  pressures  and  a  cubic  function  of  the  joint  closure 
law,  and  joint  flow  rates  are  a  linear  function  of  joint  pressure  gradients 
and  a  cubic  function  of  joint  displacements.   The  third  part  of  this  chapter 
describes  the  solution  method  (dynamic  relaxation)  used  to  solve  for  the 
nodal  unknowns. 

Two  separate  models  are  developed  in  this  thesis,  one  each  for  the 
structural  solution  and  the  fluid  solution.   However,  the  structural  solution 
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depends  on  the  fluid  solution  and  the  fluid  solution  depends  on  the 
structural  solution.   The  fourth  part  of  this  chapter  describes  the  coupling 
that  must  occur  between  the  two  solution  algorithms  for  a  truly  coupled 
solution. 

Tracer  information  can  help  immensely  when  characterizing  an 
HDR  reservoir.   Though  a  tracer  model  is  implemented,  it  was  not 
exercised  for  this  work.   The  final  section  of  this  chapter  discusses  the 
tracer  solution  implemented  by  Su  (1988). 

2.1       Element  Derivations 

This  section  describes  the  element  derivations  for  each  of  the  four 
element  types  used  in  this  thesis:  the  two  dimensional  continuum  element, 
one  dimensional  fluid  flow  element,  interface  element,  and  boundary 
spring  element.    The  two  dimensional  continuum  element  and  interface 
element  derivations  follow  those  of  Swenson  (1985);  the  derivation  of  the 
fluid  flow  element  follows  that  of  Su  (1988).   The  boundary  spring  element  is 
a  special  case  of  the  interface  element. 

2.1.1    Elasto-Dynamic  Structural  Element 

The  equation  of  equilibrium  for  a  body  in  three  dimensional 
Euclidean  space  is 

aijJ+bi=pui,  (21) 


12 


where  <Jij  is  the  stress  tensor,  b\  is  the  body  force,  p  is  the  density  and  u  is 

the  acceleration  ( indicates  the  second  derivative  of  displacement  w.r.t. 

time). 

Assuming  small  displacements,  small  strains  and  elastic  materials, 
the  strains  are  related  to  the  displacements  by 

eij=1/2(uiJ  +  uJ,i)'  (2.2) 

and  the  stresses  in  terms  of  the  strains  are 

aij  =  ^kk^ij +  2l^ij »  (2.3) 

where  5jj  is  the  Kronecker  delta,  and  X  and  u  are  Lame's  constants. 

Multiplying  the  governing  equation  2.1  by  a  small,  arbitrary  variation 
in  the  displacement,  5uj,  and  integrating  over  the  volume  we  have 

GiiiSuidV+f  bi5u;dV=  I  pu;5uidV. 

JV    JJ  Jv  jv  (2.4) 

Integrating  the  first  term  of  equation  2.4  by  parts  and  applying  the 
divergence  theorem  gives 

t:8uidS-     a::5eHdV+      biSuidV=  j  puiSiijdV, 
JS  Jv    J      J         Jv  Jv  (2.5) 

where  ti  is  the  surface  traction  over  S  and  eij  is  the  strain  tensor. 

The  finite  element  approximation  for  the  body  is  achieved  by 
discretizing  the  body  into  elements  and  applying  equation  2.5  to  each  of  the 
elements.   We  then  sum  the  contribution  of  each  element  to  obtain  the 
integrals  over  the  body: 
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I  j/s  ti5uidS-/vaIj5eijdV  +  /v  b15uidV  =  /v  pu^u^V i  . 
e=1l     e  e  )  (2.6) 

Introducing  matrix  notation  and  working  with  only  one  element  we  have 

I   5uTtdS-|    5eTadV+{    5uTbdV-|    p5uTu*dV  =  0. 
;Se  ■%  -fye  ■%  (2?) 

Over  each  element  the  displacements  are  functions  of  the  values  at  the 
nodes  surrounding  the  element 

u  =  Nua,  (2.8) 

where  ua  are  the  nodal  displacements.   The  strains  follow  as 

e  =  Bua,  (2.9) 

and  the  stresses  as 

a  =  De-  (2.10) 

Using  these  interpolation  functions  and  substituting  into  equation  2.7  yields 

I    5uaNTNdSta-       5uaBTcdV+|    5uaNTbdV 
Jr.  Jv-  Jv. 


c>q  Vg  V, 


-I    p5uaNTN  iiadV  =  0. 
^Ve  (2.11) 

Dropping  the  body  force  term,  and  since  Su^  is  arbitrary, 

I   NTNdSta-l    BTadV-l    pNTNdVu'  =0. 

JSe  a   ^e  Jv?  (212) 

Assembling  the  element  contributions,  we  obtain  the  global  stiffness 
matrix  K,  the  global  mass  matrix  M,  and  the  global  force  vector  f.  The 
nodal  accelerations  can  now  be  calculated  as: 
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where, 


u'a=  M'^f-Kua).  (213) 


Kua=    i|/vBTadV  ), 


M 


=    ij/vpNTNdv), 


ve 

(2.15) 


=  i{/sNWa) 


f  = 

(2.16) 

Note  that  Kua  is  not  evaluated  as  an  explicit  element  matrix 
multiplied  by  the  displacement  vector.  It  is  not  even  necessary  to  build  a 
global  stiffness  matrix.   If  we  place  the  external  forces  in  a  force  vector  f, 
then  accumulate  the  internal  forces  from  each  element's  contribution 
using  equation  2.13,  we  are  left  with  the  unbalanced  forces  that  result  in  the 
body  acceleration. 

2.1.2    Fluid  Flow  Element 

Figure  2.1  shows  the  differential  element  over  which  the 
conservation  of  mass  will  be  written.   The  conservation  of  mass  can  be 
stated  as, 

min-  mout=  ^(stored mass) , 
or  written  as, 
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pual 


pual  +  ^—  (pual)  dx  =  3-  (  pal)  dx  , 


3x 


dt  (2.18) 


where  p  is  the  density,  u  is  the  fluid  velocity,  a  is  the  crack  opening,  and  1  is 
the  element  thickness. 

For  a  unit  thickness  element  the  volumetric  flow  rate  is  q  =  ua.   Now 
with  constant  density  the  conservation  of  mass  can  be  written  as, 

aq    • 

5x        '  (2.19) 

Using  Darcy's  Law  we  assume  that  the  flow  rate  is  proportional  to  the 
pressure  gradient  and  the  joint  permeability, 

q  =  -k  ^ 

P8x  '  (2.20) 

The  joint  permeability,  kp,  is  given  by  the  cubic  law  (Figure  1.3), 

3 
k   =     a 

12  M  (2.21) 

where  a  is  the  joint  opening,  |i  is  the  dynamic  viscosity,  and  f  is  a  frictional 

loss  factor.  Substituting  this  into  equation  2.19  we  get  the  desired 

expression, 

d\\  Pax/        *  (2.22) 

We  now  obtain  the  weak  form  of  the  differential  equation  by 
multiplying  equation  2.22  by  an  arbitrary  weighting  function  w  and 
integrating  over  the  joint  length, 

JL    v  '  (2.23) 
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Integrating  the  first  term  by  parts  gives, 


—  ^dx 

Kp  ax  ax  ax  • 


(2.24) 
Substituting  into  equation  2.23  gives  the  desired  result, 


JL 


kn-r— -^-dx  -   |  wadx  =  0. 
P  dx  dx 


(2.25) 

The  finite  element  approximation  proceeds  as  before  by  dividing  the  joint 
into  discrete  elements.  The  integrals  of  equation  2.25  are  then  evaluated 
over  each  element  and  summed  over  the  body, 

S([*#- /*££*- twi*)-0. 

e=1V  Jl  )  (2.26) 

The  first  term  is  the  "natural"  boundary  condition  and  allows  flow  rates  to 
be  specified  where  flow  rates  are  not  equal  to  zero. 

We  introduce  shape  functions  and  matrix  notation  to  interpolate  the 
known  quantities  inside  the  element, 

a  =  N  aa, 

vt  dw     D 

W  =  Nwa,  3^=Bwa, 


ax 


q  =  N  qa, 


P  =  NPa-  a^"  =  BPa' 
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where  N  is  a  row  matrix  of  shape  functions,  B  is  the  matrix  of  the  shape 
function  derivatives  and  a,  wa,  pa,  and  qa  are  vectors  of  the  nodal  values. 
Substituting  into  equation  2.26  yields, 


/, 


BTB  padx  +   I    w  JbTB  aadx  }  =  ( 


(2.27) 
Since  wa  is  arbitrary,  and  the  nodal  values  are  constant,  we  may  write, 


i(/kpBTBdxpa  +  ( 


kDB  B  dx  pa  BTB  dx  aa)  =  0  . 


e  =  U     e  '  (2.28) 

Assembling  the  equations,  we  obtain  the  global  matrices, 

KpPa  =  Q-S  aa,  ^.29) 

where  Kp  is  the  permeability  matrix,  pa  is  the  vector  of  nodal  pressures,  Q 
is  the  vector  of  specified  flow  rates,  S  is  the  storage  matrix  due  to  the  joint 
openings  and  aa  is  the  vector  of  joint  opening  velocities  at  the  nodes.  Note 
that  the  above  formulation  recognizes  the  transient  nature  of  the  quasi- 
steady  problem  through  the  rate  of  the  joint  displacement.   Over  each  time 
interval,  the  total  flow  in  minus  the  total  flow  out  must  equal  the  change  in 
stored  fluid.   The  transients  do  not  include  inertial  terms,  but  arise  as 
changes  in  flow  rates  resulting  from  crack  opening  velocities. 

2.1.3    Interface  and  Boundary  Spring  Elements 

In  this  model,  where  we  are  approximating  rock  masses  with 
fractures  in  them,  we  need  a  way  to  model  contact  between  the  rocks  and 
the  fracture.   The  fractures  between  the  rock  masses  have  rough,  jagged 
faces  with  many  small  openings  and  pockets  that  can  contain  water  long 
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before  the  pressures  are  high  enough  to  cause  the  joint  to  open.   When  the 
fluid  pressures  are  high  enough  to  overcome  the  normal  stresses  in  the 
rock,  shear  stresses  may  cause  shearing  in  the  fracture  before  the  joint 
actually  opens. 

Swenson  (1985)  included  an  "interface  element"  that  is  implemented 
as  a  special  case  of  surface  tractions.   The  term  used  to  apply  these 
tractions  is  the  first  term  of  equation  2.12: 

I    NTNdSta, 

where  ta  are  the  tractions  at  the  nodes.  This  element  transmits 
information  through  a  joint  to  an  adjacent  element  by  using  a  specified 
traction-displacement  relation.    Figure  2.2  shows  a  typical  relationship 
which  approximates  a  rough  crack  interface  (Gangi,  1978).   Note  that  for 
zero  relative  normal  displacement  in  the  element  there  is  still  a  traction 
being  applied  across  the  crack.   This  traction  transmits  the  in-situ  stress  in 
the  rock  normal  to  the  interface  element.   A  shear  law  is  not  currently 
included,  however  provisions  are  made  so  that  a  law  such  as  that 
implemented  by  Asgian  (1988)  could  be  implemented  in  the  future.  The 
reader  is  referred  to  Swenson  (1985)  for  a  complete  description  of  the 
interface  element. 

To  include  the  far  field  stresses  in  the  model,  a  boundary  spring 
element  has  been  included.   This  boundary  spring  represents  the 
characteristics  of  the  far  field  elastic  response.   The  reservoir  is  a  local 
(although  physically  large)  perturbation  on  the  rock  surrounding  the 
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reservoir.   Using  the  boundary  springs  allows  a  small  piece  of  rock  to  be 
"cut"  out  while  retaining  the  elastic  effects  of  the  rock  outside  the  cut 
region.   This  element  is  a  simple  extension  of  the  interface  element 
discussed  above.   It  uses  a  traction-displacement  curve  to  define  a  spring 
stiffness  between  the  boundary  of  the  finite  element  mesh  and  the  far  field. 
The  boundary  spring  should  neither  be  too  soft  nor  too  stiff. 

One  approximation  assumes  that  the  local  disturbance  decays  within 
a  distance  the  same  size  as  the  reservoir.  This  leads  to  a  stress  vs 
displacement  relation  given  by, 

Au  E 

where  E  is  the  modulus  of  elasticity,  a  is  the  in-situ  stress  normal  to  the 
spring  element,  L  is  the  length  of  the  far  field  rock  mass  to  be  modeled  as  a 
spring,  and  Au  is  the  change  in  length  of  the  far  field  rock  mass  that  will 
completely  relieve  the  in-situ  stress.   Figure  2.3  shows  a  typical  stress  vs 
displacement  relationship  used  to  model  the  far  field  stresses. 

To  implement  the  finite  element  method  discussed  above,  the 
quadratic  six  noded  isoparametric  triangle  was  selected  for  the  structural 
element  (Swenson,  1985),  and  three  node  isoparametric  line  element  for  the 
fluid  and  interface  elements.    All  elements  are  integrated  numerically 
using  Gauss-Legendre  quadrature  (Zienkiewicz,  1977).   Details  of  the 
integration  for  the  structure  and  interface  elements  can  be  found  in 
Swenson  (1985).   The  fluid  element  is  similar  to  the  interface  element. 
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2.2       Joint  Model 

Two  of  the  four  element  types  described  in  this  chapter  are  used  to 
model  the  rock  joint  characteristics  as  it  is  pressurized  with  fluid.   The 
interface  element  is  superimposed  on  the  fluid  element  to  form  a 
pressurized  rock  joint.   The  surface  tractions  due  to  the  fluid  pressure  are 
added  to  the  surface  tractions  due  to  the  interface  element.   Conceptually 
this  makes  sense  because  initially,  when  there  is  no  fluid  pressure,  the 
rock  joint  carries  all  the  load.   As  the  fluid  pressure  increases,  according  to 
fundamentals  of  rock  mechanics  (Duncan,  1969),  the  load  carried  by  the 
rock  joint  decreases  until  the  fluid  pressure  is  equal  to  the  initial  stress  in 
the  rock.   Using  rock  mechanics  concepts  this  is  stated  as: 

Total  stress    =    Effective  stress    +    Pore  pressure 
(Joint  stress)  (Fluid  pressure) 

The  pore  pressure  is  simply  the  fluid  pressure  and  is  applied  as  a  surface 

traction.  The  effective  stress  is  specified  by  a  stress  vs  displacement  curve 

that  is  called  the  opening  law  and  is  also  applied  as  a  surface  traction.   The 

opening  law  used  in  this  thesis  is  the  "Bed-of-Nails"  model  (Gangi,  1978). 

Gangi  showed  that  the  functional  dependence  of  the  joint  opening  variation 

of  a  fracture  can  be  modeled  as 

m 
1  -" 


a(P)  =  a0 

°L      \/rV    J'  (2.30) 

where  a0  is  the  zero  pressure  joint  opening,  P  is  pore  pressure,  Pi  is  the 
effective  modulus  of  the  asperities,  and  m  is  a  constant  (0  <  m  <  1)  which 
characterizes  the  distribution  function  to  the  asperity  lengths.   For  all  work 
in  this  thesis  a0  =  0.32  mm,  Pi  =  70  MPa  and  m  =  0.3636  (Brown,  1988a). 
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Figure  2.4  compares  the  "Bed-of-Nails"  model  with  a  "natural" 
crack.    Figure  2.4a  shows  a  "natural"  crack  formed  by  creating  a  hairline 
fracture,  then  translating  the  lower  half  of  the  medium  to  the  right.   This  is 
similar  to  how  open  cracks  are  formed  in  nature.   Figure  2.4b  illustrates 
the  "Bed-of-Nails"  model.   The  distribution  of  asperities  is  treated  as  a 
distribution  of  rods,  which  is  much  simpler  to  analyze  than  the  mechanical 
properties  of  the  natural  crack. 

The  joint  law  used  in  the  present  calculations  is  very  stiff.   As  a 
result,  it  was  necessary  to  use  care  to  ensure  convergence.   This  is 
discussed  in  Section  3.1  and  Chapter  VI. 

2.3       Solution  Method 

As  mentioned  before,  the  problem  we  are  solving  is  highly  nonlinear. 
Many  methods  exist  for  solving  highly  nonlinear  problems,  but  few  are  as 
robust  as  Dynamic  Relaxation  (DR).   Consider  a  simple  problem  where  a 
stiffness  matrix  is  a  function  of  displacement.   The  static  equilibrium 
equation  can  be  written  as, 

K(x)x  =  f,  (2.31) 

where  K(x)  is  the  nonlinear  stiffness  matrix,  x  is  the  displacement  vector, 
and  f  is  the  force  vector.  To  solve  this  problem  using  DR  we  would  rewrite 
the  equation  as, 

Mx  +  Cx  +  Kx  =  f,  (2.32) 
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where  M  is  a  mass  matrix  and  C  is  a  damping  matrix.    If  the 
damping  matrix  is  defined  as  C  =  cM,  then  equation  2.32  can  be  rewritten 
as 

x=M"  (f  -  Kxj-cx.  (2.33) 

When  the  mass  matrix  is  lumped  to  be  a  matrix  with  elements  on  the 
main  diagonal,  the  inverse  of  the  mass  matrix  is  the  inverse  of  each  of  the 
nodal  masses.   As  discussed  in  Section  2.1.1,  an  equation  of  the  form  of 
equation  2.33  is  vectorizable  and  does  not  require  that  the  entire  stiffness 
matrix  be  assembled.   These  features  greatly  reduce  storage  requirements 
and  simplify  programming  requirements.    Note  also  that  the  damping 
term  is  a  scalar,  not  a  vector  or  matrix  quantity.  Equation  2.33  can  now  be 
integrated  explicitly  by  the  central  difference  method.   The  resulting 
equations  are: 

xn  +  i/  =  x_i/  +  xnAt , 

n+/2         n    /2         n  (2M) 


n+1      xn-l 


Xn+l/At 


'2  (2.35) 

where  n  is  the  time  step  and  At  is  a  fixed  pseudo-time  increment.   Now  we 
select  At,  M,  and  c  to  allow  the  solution  to  converge  as  fast  as  possible. 

The  solution  algorithm  for  the  structure  relies  heavily  on  the  original 
algorithm  in  Swenson's  (1985)  code.   Because  it  is  a  dynamic  code,  the 
inertia  and  physical  mass  matrix  already  exist.    We  alter  this  original 
mass  matrix  to  maintain  faster  convergence.    The  density  for  each  element 
is  divided  by  the  square  of  the  information  transit  time  for  that  element. 
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This  has  the  effect  of  setting  the  transit  time  to  ~1  for  each  element,  which 
also  sets  the  minimum  integration  time  step  to  ~1.  Underwood  (1983) 
shows  an  example  problem  that  describes  this. 

For  the  fluid  solution  an  explicit  stiffness  matrix  is  built  on  the 
element  level.   When  this  is  available,  Underwood  (1983)  suggests  using 
Gerschgorin's  theorem  to  make  the  mass  matrix  a  function  of  the  stiffness 
matrix: 

mii  =  1/4  At2  Zj  I  kij  I  ,  (2.36) 

where  m^  are  the  diagonal  elements  of  the  diagonal  pseudo  mass  matrix 
M,  and  kjj  are  the  elements  of  the  stiffness  matrix  K.   He  also  suggests 
evaluating  M  with  At  =  1.1  and  iterating  with  At  =  1.0  to  ensure  stability.  Su 
(1988)  found  it  necessary  to  iterate  at  At  =  0.5  for  stability. 

Underwood  (1983)  also  suggests  using  Rayleigh's  quotient  to  predict 
the  approximate  minimum  natural  frequency  (0o  for  the  current 
deformation  mode, 


co0  =  V  xTR  x  /  xTM  x  .  (2.37) 

Then  the  damping  for  this  mode  is  approximated  by 

c  =  2^Wo .  (2.38) 

where  t,  is  the  damping  ratio. 

The   damping  ratio  %  can  be  used  to  reduce  or  increase  the  damping 
factor  for  better  convergence.  Because  0)o  is  only  an  approximation  to  the 
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lowest  active  frequency  there  is  only  a  general  correlation  between  the 
damping  ratio  and  the  rate  of  convergence. 

This  method  works  reasonably  well  for  both  the  pressure   and 
structure  solutions  as  long  as  no  rigid  body  modes  are  present  in  the 
solution.   Rigid  body  displacements  correspond  to  low  frequencies,  and  the 
damping  calculated  by  equation  2.38  becomes  extremely  small.   As  a  result, 
higher  frequencies  are  underdamped  and  convergence  is  very  slow,  with 
much  high  frequency  oscillation.    For  the  problems  worked  in  this  thesis 
rigid  body  modes  are  almost  always  present  in  both  the  structure  and  the 
fluid.   Therefore,  an  alternate,  less  sophisticated,  but  more  reliable  velocity 
reduction  damping  method  was  implemented.    After  the  accelerations  are 
integrated  to  update  the  velocities,  the  new  velocities  are  multiplied  by  a 
reduction  or  damping  term, 


*n  +  y2=  Un-l/Z  +  ^AtJ  damp 


(2.39) 

This  slowly  reduces  the  velocity  as  the  problem  progresses  and  eventually 
forces  the  problem  to  steady  state. 

The  major  problem  with  this  approach  is  that  the  appropriate 
damping  is  problem  dependent.   This  places  a  greater  burden  on  the  user  to 
examine  results  closely  to  ensure  convergence.    However,  due  the 
interactive  nature  of  the  program,  the  analysis  can  be  stopped  at  any  time 
and  the  response  monitored  to  allow  judgment  to  dictate  whether  the 
damping  should  be  increased  or  decreased.   This  allows  the  user  to  supply 
information  about  how  close  the  solution  is  to  steady  state  that  would 
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automatically  be  supplied  by  a  method  such  as  Raleigh's  quotient  if  the 
rigid  body  modes  were  not  present. 

The  value  of  damping  should  not  be  taken  lightly  however.   From 
equation  2.39  we  can  see  that  as  xn+i/2  is  substituted  for  xn_i/2  at  each  time 
step  the  velocities  are  reduced  as  a  function  of  Cn  where  n  is  the  number  of 
iterations.  For  n  =  2000  the  effect  of  damping  on  the  velocity  can  readily  be 
seen: 

Damping       Velocity  Reduction 


0.999 

0.135 

0.9999 

0.819 

0.99999 

0.980 

From  this  we  can  see  that  small  changes  in  the  damping  parameter 
can  have  profound  effects  on  the  solution  over  several  thousand  iterations. 
If  the  damping  factor  is  increased  slightly  the  velocity  reduction  will  not 
occur  so  readily  and  the  solution  can  progress  to  convergence  at  a  faster 
rate.   It  is  the  user's  responsibility  to  recognize  these  inconsistencies  while 
analyzing  results  and  to  modify  parameters  such  as  the  damping  factor 
and  tolerance  and  iteration  coupling  parameters  (discussed  in  Section  3.1) 
to  ensure  that  the  results  obtained  are  understood  thoroughly. 

2.4       Coupling  between  Fluid  and  Structural  Models 

As  shown  above,  both  the  fluid  and  structural  equations  are 
independent  modules.   The  structural  model  is  coupled  to  the  fluid  model 
through  the  pressures  applied  on  the  joints.   The  fluid  model  is  coupled  to 
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the  structural  model  through  the  joint  opening  displacements.    In  the  fluid 
solution,  complete  compatibility  between  the  models  is  obtained.   Both  the 
models  are  solved  iteratively.   Information  is  passed  between  the  models  at 
user  specified  intervals.    This  process  continues  until  global  convergence  is 
obtained. 

2.5       Tracer  Model 

Tracers  are  used  to  track  the  motion  of  fluid  inside  an  HDR  reservoir. 
Tracer  studies  can  help  determine  reservoir  volume  and  thermal 
characteristics  as  well  as  other  information.   At  Los  Alamos  both  reactive 
and  nonreactive  tracers  have  been  used.   Su  (1988)  developed  an  uncoupled 
tracer  model  for  nonreactive  tracers.   His  tracer  model  has  not  been 
exercised  in  this  work,  however  it  is  included  in  the  model.   Some 
modifications  may  be  needed  before  it  is  used  extensively.  Robinson  (1988) 
has  shown  that  the  joint  aperture  appropriate  for  the  fluid  flow  solution 
may  not  be  appropriate  for  the  tracer  fluid  flow  solution.   Figure  2.5 
compares  the  hydraulic  aperture  wh  with  the  tracer  aperture  wt.  The 
tracer  aperture  is  much  larger  than  the  hydraulic  aperture  used  to  find  the 
fluid  flow  solution.   The  larger  aperture  for  the  tracer  leads  to  larger 
volumes  and  significantly  influences  tracer  concentrations.    Robinson 
(1985)  also  developed  a  residence  time  distribution  curve  which  is  useful  for 
determining  a  reservoir's  volume.   His  method  could  be  used  to  verify  the 
tracer  model. 
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Figure  2.1:       Conservation  of  Mass  Differential  Element 
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Figure  2.2:       Typical  Joint  Opening  Law 

Traction  —  Displacement  Relation 
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Figure  2.3:       Typical  Boundary  Spring 

Traction  —  Displacement   Relation 
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A)  Natural  Crack   (stylized) 
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B)  Mechanically  and  Hydraul  icaily  Equivalent  Crack 
(schematic ) 


Figure  2.4:       "Bed  of  Nails"  Joint  Model 
(From  Gangi,  1978) 
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Figure  2.5:       Hydraulic  Aperture  and  Tracer  Aperture 
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Chapter  III 
Implementation  and  Modelling  Approach 

A  computer  program  called  DRACULA  (Dry  £ock  Analysis  Code) 
has  been  developed  to  simulate  transient  fluid  flow  in  fractured, 
deformable,  rock  masses.    This  code  implements  the  equations  and 
concepts  discussed  in  Chapter  II.   This  chapter  discusses  the  basic 
concepts  necessary  for  operation  of  DRACULA,  then  presents  some  general 
information  about  modelling  HDR  reservoirs. 

3.1       Implementation 

DRACULA  was  developed  to  provide  an  interactive  graphics 
environment  to  allow  the  user  to  specify  and  monitor  the  problem  in  a 
simple,  reliable  way.   To  monitor  the  solution's  progress,   the  user  can  stop 
an  analysis  and  view  intermediate  results.    If  necessary,  control 
parameters  can  be  modified  before  restarting  the  analysis.   The  code 
CRACKER  (Swenson,  1985)  was  used  as  a  basic  foundation  to  implement 
these  interactive  concepts  and  to  make  use  of  the  existing  structural 
analysis  concepts.   The  interactive  nature  of  CRACKER  unifies  the 
traditionally  separate  tasks  of  preprocessing,  analysis,  and  post- 
processing. 

From  the  MAIN  menu  page  the  user  may  define  the  problem,  save 
data,  run  the  analysis  or  view  the  results  (see  Fig  3.1).   The  ANALY 
PARAM  option  controls  three  options:  the  GLOBAL,  FLUID  and  DYN  REL 
parameter  pages.    Each  of  these  contains  parameters  that  control  program 
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execution  such  as  end  time,  time  step  size,  data  output  frequency  and 
analysis  type.  The  PLOT  page  allows  the  user  to  view  the  results  as 
pressure  and  stress  contour  plots,  time  history  plots,  displaced  mesh  plots 
and  others.   SAVE  RSTRT,  SAVE  PLOT  and  SAVE  GEO  allow  the  user  to 
output  the  restart,  plot  and  geometry  data  to  files  for  future  reference.  The 
REMESH  page  allows  the  user  to  delete  and  add  elements,  drag  nodes, 
change  material  properties,  and  redefine  initial  and  boundary  conditions. 
GO  BATCH  and  GO  INTERACT  tell  the  program  to  start  a  batch  analysis 
or  perform  an  interactive  analysis.   For  the  batch  analysis,  a  batch  start  file 
is  written  and  the  program  stops,  ready  to  restart  in  batch  mode.  For  an 
interactive  analysis,  as  the  analysis  proceeds  messages  are  displayed  to  the 
user  in  the  MAIN  PAGE.   If  trouble  is  encountered  a  message  is  displayed 
and  control  returns  to  the  user;  in  batch  mode  if  trouble  is  encountered  a 
restart  file  is  written. 

Three  of  the  control  parameters  contained  in  the  DYN  REL  and 
FLUID  parameter  pages  regulate  the  total  number  of  iterations  the  problem 
will  be  allowed  to  run  as  well  as  how  often  coupling  occurs  between  the  two 
solutions.    In  the  FLUID  page  the  "Number  of  Fluid  coupling  iter"  is  the 
number  of  iterations  the  fluid  routine  is  allowed  before  it  returns  the 
current  solution  to  the  structure.   In  the  DYN  REL  page  the  "Number  of 
Structure  coupling  iter"  is  the  number  of  iterations  the  structure  routine  is 
allowed  before  it  returns  the  current  solution  to  the  fluid.   In  the  DYN  REL 
page  the  "Total  allowed  iterations  (structure)"  is  the  total  allowed  iterations 
for  the  problem.   It  is  designated  for  the  structure  because  it  actually  counts 
the  total  iterations  for  the  structure  solution.  This  is  because  the  structure 
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iterations  are  more  time  consuming  than  fluid  iterations.    Also,  the  fluid 
solution  typically  converges  in  less  iterations  than  the  structure  and  can 
therefore  usually  have  a  better  fluid  solution  for  every  set  of  structure 
iterations. 

Two  other  control  parameters  contained  (one  each)  in  the  FLUID  and 
DYN  REL  pages  are  the  convergence  tolerances  for  each  of  the  two 
routines.   The  fluid  routine  uses  a  different  convergence  criterion  than  the 
structure  routine.    The  fluid  routine  is  converged  when  each  individual 
normalized  pressure  change  from  the  last  iteration  to  this  iteration 
changed  less  than  its  tolerance: 

"new"  "old     .  „  , 

5 <  to1- 

pold 

The  structure  is  converged  when  the  square  root  of  the  normalized  sum  of 
the  displacements  squared  is  less  than  its  tolerance: 
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In  a  reservoir  there  are  many  fractures  that  are  able  to  store  fluid 
before  building  enough  pressure  to  relieve  stresses  holding  the  cracks 
together.   To  model  this  phenomenon  the  fluid  elements  have  a  special 
characteristic  —  they  are  not  required  to  start  with  fluid  in  them.   If  the 
fluid  element's  initial  opening  (a  fluid  material  property)  is  greater  than  a 
user  specified  tolerance,  then  the  element  is  assumed  to  be  empty.   When 
adjacent  fluid  elements  fill  with  fluid  they  can  start  to  fill  the  next  empty 
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element,  and  the  fluid  solution  can  progress  through  the  mesh.   It  is  not 
necessary  that  any  or  all  of  the  fluid  elements  be  open,  or  closed  initially. 
By  closing  them  all,  the  steady  state  flow  solution  for  the  reservoir  can  be 
reached  at  the  end  of  the  first  converged  solution.  By  opening  them  all,  the 
fluid  movement  through  the  reservoir  can  be  observed  from  the  time  fluid 
injection  begins. 

The  dynamic  relaxation  solution  method  is  implemented  for  both  the 
fluid  and  structural  finite  element  schemes  in  separate  subroutines.   The 
structure  routine  solves  the  structural  finite  element  problem,  controls 
structure  data  output  to  files,  and  controls  the  coupling  between  itself  and 
the  fluid  routine.  The  fluid  routine  solves  the  fluid  finite  element  problem, 
controls  fluid  data  output  to  files,  and  presents  the  structure  routine  with 
an  updated  fluid  solution. 

The  time  step  reduction  factor  is  another  important  parameter 
contained  in  the  GLOBAL  parameter  page.   The  solution  method  currently 
implemented  is  an  explicit  integration  method.   The  stability  of  the 
integration  scheme  is  dependent  on  the  time  step  size.   The  solution 
algorithm  will  select  an  appropriate  time  step  based  on  the  size  of  the 
smallest  structural  element,  but  the  stiffnesses  of  the  interface  elements 
and  specifically  the  joint  law  also  need  to  be  considered.   Because  the  joint 
law  is  very  stiff  compared  to  a  typical  structural  element  for  these 
problems,  the  time  step  must  be  reduced  by  a  factor  of  100.  This  increases 
the  run  time  on  a  problem  because  100  times  as  many  time  steps  are  now 
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required  for  the  same  displaced  solution.    Chapter  VI  gives 
recommendations  on  removing  this  constraint. 

3.2       Modelling  Approach 

HDR  reservoirs  are  large,  three  dimensional  underground  regions, 
typically  measured  in  hundreds  of  meters.    These  three  dimensional 
regions  comprise  a  lattice  of  small  and  large  fractures.    When  fluid  is 
pumped  into  a  reservoir  the  reservoir  expands  like  a  balloon.   Most  of  the 
fluid  can  be  recovered  by  letting  the  reservoir  contract,  but  some  of  the  fluid 
will  be  lost  at  the  outer  boundary  of  the  reservoir  and  some  is  trapped  in 
"open"  fractures.   The  Phase  II  HDR  reservoir  at  Fenton  Hill,  New  Mexico 
displays  these  characteristics.  As  a  first  approximation  the  Phase  II 
reservoir  will  be  modelled  as  a  horizontal  plane  with  unit  depth.  We  look  at 
the  horizontal  plane  in  plan  view,  with  flow  occurring  only  in  the  fractures 
in  this  plane. 

The  researchers  at  Los  Alamos  National  Laboratory  have  evidence 
that  as  99%  or  more  of  the  volume  in  their  HDR  reservoir  is  due  to  the 
small,  low  opening  stress  fractures  (Brown,  1988).   A  small  number  of  the 
fractures  in  the  reservoir  have  a  high  opening  stress.   The  low  opening 
stress  fractures  are  called  "tensile"  fractures,  and  the  high  opening  stress 
fractures  are  called  "shear"  fractures.   Brown  (1988)  also  reports  that  the 
principal  stresses  on  the  reservoir  are  about  1 0  MPa  and  24  MPa  for  the 
minimum  and  maximum  principal  in-situ  stresses.    For  this  model  the 
shear  fractures  are  assumed  to  be  perpendicular  to  the  maximum 
principal  stress  and  the  tensile  fractures  are  assumed  to  be  perpendicular 
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to  the  minimum  principal  stress.    In  reality  the  principal  stresses  would  be 
rotated  counterclockwise  slightly,  but  DRACULA  does  not  currently  have  a 
shear  law  in  the  joint  model.   The  in-situ  stress  on  the  tensile  fractures  is 
10  MPa,  and  the  in-situ  stress  on  the  shear  fractures  is  24  MPa;  both  are 
strictly  normal  stresses. 

3.2.1    Joint  Model  Implementation 

As  discussed  in  Section  2.2,  the  joint  model  uses  the  fluid  element 
superimposed  on  the  interface  element.   The  tensile  fractures  and  shear 
fractures  use  the  same  joint  model  but  refer  to  different  materials.   The 
tensile  fractures  use  material  type  one;  the  shear  fractures  use  material 
type  two.   The  fluid  and  interface  elements  get  their  material  data  from 
their  respective  material  tables.   Figure  3.2  shows  the  joint  opening  law  for 
interface  types  one  and  two.   There  is  only  one  opening  law  as  defined  in 
Section  2.2.  The  curve  is  shifted  to  keep  the  body  in  initial  equilibrium  for 
the  two  in-situ  stresses.  The  distance  the  curve  is  shifted  must  be  stored  as 
fluid  material  data  to  be  used  as  an  initial  opening  for  the  joints.  This  gives 
each  joint  some  finite  free  volume  before  any  fluid  is  pumped  in. 

Figure  3.3  shows  the  finite  element  mesh  used  for  the  problems 
analyzed  in  Chapter  IV.   The  three  vertical  joints  model  the  shear 
fractures.   The  tensile  joints  are  offset  to  force  the  flow  paths  to  include  the 
high  opening  stress  shear  joints. 
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3.2.2    Reservoir  Boundary  Conditions 

To  model  the  effects  of  the  rock  surrounding  the  reservoir,  springs 
are  placed  on  the  outer  edge  of  the  finite  element  mesh.   Figure  3.4  shows 
typical  plots  for  interface  materials  three  and  four.   Notice  that  the  springs 
are  preloaded  to  10  MPa  and  24  MPa  to  maintain  initial  equilibrium.   As  the 
reservoir  expands  and  the  volume  increases  the  far  field  stress  increases 
with  the  spring's  compression.    The  spring  element  sets  the  displacements 
on  the  outside  of  the  element  to  zero  and  has  inside  nodes  attached  to  the 
rock  masses. 

The  fluid  boundary  conditions  are  important  to  consider  before 
running  a  problem.   Pressure  and  flow  rate  histories  may  be  specified  for 
transient  problems.    The  "natural"  boundary  condition  for  the  fluid  flow 
problem  is  that  no  flow  occurs  on  a  boundary  where  no  condition  is 
specified.  One  other  condition  that  can  be  applied  globally  to  the  reservoir 
boundary  is  a  condition  that  allows  fluid  to  "leak"  out  of  the  reservoir  as  a 
function  of  flow  rate  vs.  pressure.   This  is  called  a  far  field  flow  loss  or 
leakage  boundary  condition. 


39 


MAIN 
PAGE 


ANALY  PARAM 


PLOTS 


SAVE  RSTRT 


SAVE  PLOT 


SAVE  GEO 


REMESH 


GO  BATCH 


GO  INTERACT 


REMESH 
PAGE 

DRAG 

NODE 

DEL 

ELMT 

ADD 

ELMT 

BOUND 

COND 

MAT 

PROP 

INIT 

COND 

PRINC 

VECT 

ANALY 
PARAM 


GLOBAL 


FLUID 
DYN  REL 


FINISHED 


PLOT 
PAGE 


DISPLACE 


VELOCITY 


PRINC  VECT 


SIG  3D  CNTR 


SIG  2D  CNTR 


EPS  CONTR 


LINE  PLOT 


TIME  HISTRY 


FLUID 


BOUND 
COND 

DISPL 

BC 

VEL 

BC 

PRESS 

NORM 

PRESS 

SHR 

POINT 

LOAD 

PRESS 

vs  t 

POINT 

vs  t 

FLUID 

BC 

Figure  3.1:       DRACULA  Menu  Concept 
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Figure  3.2:       Joint  Opening  Law, 
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Figure  3.3:       Application  Problems  Finite  Element  Mesh 
(Note:  Only  flow  paths  shown  for  clarity) 
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Chapter  IV 
Verification  Problems 

Two  verification  problems  will  be  used  to  illustrate  DRACULA's 
ability  to  model  fluid  flow.  A  single  flow  path  is  used  so  the  results  can  be 
interpreted  more  easily.   The  path  considered  is  one  tensile  flow  path  from 
Figure  3.3.  The  first  problem  is  a  steady  state  flow  problem.  A  specified 
flow  rate  is  applied  at  one  boundary  of  the  problem  and  a  far  field  flow 
boundary  condition  is  applied  to  allow  fluid  leakage  at  the  other  boundary. 
The  second  problem  is  a  transient  flow  problem.   A  constant  reservoir 
pressure  is  applied  and  the  steady  state  (no  fluid  flow)  solution  is  found. 
Then  a  step  change  in  the  pressure  boundary  condition  initiates  a  pressure 
transient. 

Figure  4.1  shows  the  mesh  used  for  both  of  the  problems  discussed  in 
this  chapter.    The  verification  problems  are  simple,  one  dimensional 
problems  with  a  single  flow  path  in  the  X  direction.   Table  4.1  summarizes 
the  reservoir  properties  used  for  all  the  problems  in  this  thesis.   Poisson's 
ratio  was  set  to  zero  for  the  verification  problems.  Figure  4.2  shows  plots  of 
the  boundary  spring  stiffness  applied  for  both  verification  problems. 
Boundary  springs  with  a  relatively  high  stiffness  were  used  to  help  speed 
convergence. 

4.1       Steady  State  Verification  Problem 

This  problem  displays  the  option  to  specify  a  pressure  -  flow  rate 
relation  on  the  boundary  of  the  mesh  to  simulate  a  far  field  "leakage" 


44 


boundary  condition.   Figure  4.3a  shows  the  far  field  relation  specified  as  a 
linear  function  of  pressure  and  flow  rate  (nonlinear  functions  can  also  be 
specified).  The  flow  rate  was  specified  at  the  right  end  of  the  mesh  as 
shown  by  Figure  4.3b.   Normally  the  far  field  pressure  -  flow  rate  relation 
would  be  enforced  on  the  right  boundary  as  well  as  the  left  boundary.  But 
when  a  specific  pressure  or  flow  rate  condition  is  specified  on  a  boundary, 
the  far  field  boundary  condition  will  not  be  applied  at  that  node. 

The  flow  rate  solution  should  be  constant  through  the  entire  flow 
path.   Also  note  that  numerically  all  the  boundary  conditions  specified  for 
this  problem  are  flow  rates;  however,  because  the  far  field  boundary 
condition  is  a  function  of  pressure,  a  unique  pressure  solution  does  exist. 
Convergence  rates  are  normally  slower  for  problems  with  far  field  and  flow 
rate  boundary  conditions  specified. 

Figures  4.4  is  a  flow  rate  plot  for  this  steady  state  problem.  The  flow 
rate  is  constant  through  the  flow  path  as  expected.   As  shown  by  Figure  4.5, 
the  pressure  on  the  left  boundary  corresponds  to  the  correct  flow  rate  for  the 
far  field  boundary  condition.   In  addition,  the  joint  stress  shown  by  Figure 
4.6  added  to  the  fluid  pressure  from  Figure  4.5  gives  the  total  stress  shown 
by  the  line  plot  of  the  Y  stress  in  Figure  4.7.  This  shows  that  the 
fundamental  rock  mechanics  assumption  from  Section  2.2  is  satisfied  in 
DRACULA. 

The  last  figure  for  this  problem  (Figure  4.8)  is  a  displaced  mesh  plot. 
Note  that  the  joint  openings  change  as  the  pressure  changes  in  the  joint. 
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Since  the  flow  rate  is  constant  through  the  flow  path,  the  pressure  gradient 
must  also  adjust  with  the  joint  openings  to  maintain  the  specified  flow  rate. 

4.2       Transient  Verification  Problem 

This  problem  is  analogous  to  inflating  a  balloon.   All  the  input  flow  is 
used  to  expand  the  joint.   If  the  pressure  is  reduced,  this  stored  fluid  is 
available  to  flow  out  of  the  joint.  To  start  the  problem  a  constant  9.0  MPa 
fluid  pressure  was  applied  to  one  end  of  the  flow  path  and  a  zero  flow 
boundary  condition  was  imposed  on  the  opposite  end  of  the  flow  path.  This 
results  in  a  fluid  solution  with  constant  fluid  pressure  in  the  joint  and  no 
fluid  flow.   Next  the  applied  pressure  is  dropped  to  1 .0  MPa,  initiating  a 
pressure  transient  in  the  problem,  causing  fluid  flow  out  of  the  joint. 

The  analytical  solution  for  the  first  part  of  this  problem  comes  from  a 
simple  one  dimensional  rock  compression  problem.    The  fluid  pressure  is 
constant  in  the  joint,  causing  a  uniform  linear  translation  of  the  two  rock 
masses.  The  initial  stress  in  the  block  is  ax  =  -24  MPa  and  oy  =  -10  MPa. 
Initially  the  joint  also  has  a  stress  of  -10  MPa.  Inflating  the  joint  to  a  fluid 
pressure  of  9  MPa  opens  the  joint.  This  causes  the  block  to  translate  and 
increases  the  boundary  load  to  -10.2  MPa.   In  the  joint,  the  fluid  pressure  is 
9  MPa  and  the  joint  stress  is  reduced  to  1.3  MPa.  The  stress  in  the  block 
also  increases  to  10.3  MPa.  The  values  are  not  exact  because  relatively  loose 
convergence  tolerances  were  used.   The  effect  of  this  can  be  seen  on  the  line 
plot  in  Figure  4.9.   The  Y  stress  should  be  constant  along  this  line  plot. 
Figure  4.10  is  a  time  history  plot  that  also  shows  that  though  the  solution  is 
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nearly  converged,  more  iterations  are  needed  to  allow  the  blocks  to  translate 
and  increase  the  boundary  spring  load  to  -10.3  MPa. 

At  time  =  0+  the  applied  pressure  was  dropped  to  1.0  MPa.    Each 
transient  solution  must  be  obtained  by  iterating  on  the  nonlinear  flow 
problem,  with  the  previous  solution  giving  initial  conditions  for  the  solution 
currently  being  sought.   An  important  parameter  that  must  be  selected 
before  a  transient  is  initiated  is  the  time  step  (t-step)  over  which  the 
transient  will  act.   Recall  from  Chapter  II  that  the  fluid  flow  is  governed  by, 

KpPa  =  Q-S  aa. 

The  S  aa  term  acts  as  a  fluid  source  in  the  analysis.   If  a  joint  is  closing, 
this  term  supplies  additional  flow  to  the  calculation.   In  DRACULA  the  aa 
term  is  implemented  in  finite  difference  form  as, 

a  new "  a  old 


aa= 


t-step 


The  fluid  flow  equations  we  are  solving  are  now  a  linear  function  of 
the  joint  displacements  from  the  mesh  (anew)  and  still  a  cubic  function  of 
joint  displacements  in  the  joint  permeability.   The  flow  rates  are  influenced 
directly  by  the  selection  of  t-step  as  well  as  the  displacements.  Notice  above 
that  as  t-step  increases,  the  storage  term  has  less  and  less  effect  on  the 
solution,  and  approaches  the  steady  state  solution.   And  conversely,  as  t- 
step  is  reduced,  the  flow  rates  can  be  influenced  significantly  by  the  storage 
term. 
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Figure  4.11  shows  a  time  history  of  four  displaced  mesh  plots  for  a 
transient  solution  with  t-step  =  0.125.  As  shown  in  Figure  4.12,  the  initial 
condition  at  time  =  0.0  is  a  uniform  joint  opening  with  a  pressure  of  9.0 
MPa.  When  the  pressure  at  the  outlet  is  dropped,  fluid  begins  to  flow  out  of 
the  joint.  The  joint  closes  first  at  the  outlet.  Eventually,  the  excess  fluid  is 
forced  out  of  the  joint  and  the  pressure  is  uniform  at  1.0  MPa. 

Figure  4.14  is  a  history  plot  of  the  fluid  flow  rate  at  the  outlet  and  a 
history  plot  of  the  pressure  at  the  opposite  end  of  the  flow  path.  The  outlet 
flow  rate  rapidly  decreases  as  the  pressure  driving  the  flow  decreases  to  the 
specified  outlet  pressure. 

Figure  4.15  is  a  similar  fluid  flow  rate  plot  comparing  the  outlet  flow 
rate  for  different  values  of  t-step.  This  plot  displays  the  effects  t-step  has  on 
the  storage  term  as  discussed  above.   Smaller  values  of  t-step  cause  higher 
initial  flow  rates  with  sharper  transients;  larger  values  of  t-step  cause 
lower  initial  flow  rates  with  solutions  nearer  a  new  steady  state  solution  for 
the  new  boundary  condition. 

One  simple  check  for  the  validity  of  a  transient  solution  for  this  type 
of  problem  is  that  the  integral  over  time  of  the  flow  out  of  the  joint  be  equal  to 
the  change  in  volume  of  the  joint.   Table  4.2  summarizes  the  changes  in 
volume  representing  the  integrals  for  the  area  under  each  flow  rate  vs. 
time  curve  from  Figure  4.15.   As  t-step  decreases  from  0.5,  the  volume 
integrated  under  the  flow  rate  curve  approaches  the  actual  change  in 
volume.  Note  that  for  t-step  =  1.0  the  integrated  volume  is  also  increasing 
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and  approaching  the  actual  change  in  volume.   As  the  time  step  is 
increased  the  solution  does  approach  the  next  steady  state  solution. 

4.3       Summary 

In  this  chapter  we  presented  two  verification  problems  to  illustrate 
DRACULA's  ability  to  model  coupled  fluid  flow  and  rock  deformation  in  an 
HDR  reservoir.   The  first  problem  displayed  steady  state  results  from  flow 
rate  and  far  field  fluid  boundary  conditions.   The  finite  element  results 
show  that  the  flow  rate  is  constant  through  the  flow  path  at  the  correct 
specified  value  and  that  the  pressure  solution  does  satisfy  the  far  field 
boundary  condition.  The  second  problem  showed  the  effects  of  the  time  step 
size  on  the  results  of  a  transient  solution.   This  problem  shows  that  the  time 
step  size  does  affect  the  results  of  a  transient  problem  and  that  as  the  time 
step  size  is  reduced  a  more  accurate  transient  solution  is  obtained.   Both 
problems  show  that  the  effective  stress  law  from  Section  2.2  is  satisfied. 
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Table  4.1 
Reservoir  Properties 


Poisson's  Ratio  0.20 

Young's  Modulus  25000  MPa 

Initial  Joint  Aperature,  x  0.1032e-3  m 

Initial  Joint  Aperature,  y  0.1632e-3  m 

Maximum  Principal  Stress,  ax  24  MPa 

Minimum  Principal  Stress,  ay  10  MPa 

Fluid  Dynamic  Viscosity,  u  (@  220°  C)  11 6.6e-6  N-s/m2 

Factor  of  Roughness,  f  1.5 
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Table  4.2 


Joint  Volume  Change 


t-step  (s) 

0.0625 

0.125 

0.25 

0.5 

1.0 

oo 

Vol.  (m3x106) 

933.5 

925.9 

916.5 

914.7 

947.3 

957.9 
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Structural  Element 


Fluid  Element 


Interface  Element 


Figure  4.1:       Verification  Problems  Finite  Element  Mesh 
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Figure  4.2:       Verification  Problems  Boundary  Springs 
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Figure  4.3a:     Far  Field  Flow  Loss  Boundary  Condition 
(Steady  state  verification  problem) 
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Figure  4.3b: 


Flow  Rate  Boundary  Condition 
(Steady  state  verification  problem) 
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Figure  4.4:       Flow  Rate  Contour  Plot 

(Steady  state  verification  problem) 
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Figure  4.5:       Pressure  Contour  Plot 

(Steady  state  verification  problem) 


Figure  4.6:       Joint  Effective  Opening  Stress 

(Steady  state  verification  problem) 
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Figure  4.7a:     Y  Stress  Line  Plot 

(Steady  state  verification  problem) 
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Figure  4.7b:     Line  of  Application  for  Above  Line  Plot 
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Figure  4.8:       Displaced  Mesh  Plot 

(Steady  state  verification  problem) 
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Figure  4.9a:     Y  Stress  Line  Plot,  time  =  0.0 

(transient  verification  problem) 
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Figure  4.9b:     Line  of  Application  for  Above  Line  Plot 


59 


NODE  NUMBER   28 


-0.600-"- 

0.000   0.500   1.000   1.500   2.000   2.500 


TIME 


(xlO 


Figure  4.10:     Time  History  Plot,  Initial  9  MPa  Solution 
(transient  verification  problem) 
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Figure  4.11: 


Time  History  Plot  of  Displaced  Mesh  Plots, 
t-step  =  0.125  (transient  verification  problem) 
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Figure  4.11  cont: 


Time  History  Plot  of  Displaced  Mesh  Plots, 
t-step  =  0.125  (transient  verification  problem) 
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Figure  4.12: 


Initial  Displaced  Mesh  Plot 
(transient  verification  problem) 
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Figure  4.13:     This  figure  intentionally  left  blank 
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Figure  4.14a:  Flow  Rate  Transient  Solution 

(transient  verification  problem) 
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Figure  4.14b:   Pressure  Transient  Solution 

(transient  verification  problem) 
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Figure  4.15:     Fluid  Flow  Rate  Solution  for  Different 

values  of  t-step  (transient  verification  problem) 
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Chapter  V 
Applications 

This  chapter  presents  the  results  of  several  problems  analyzed  to 
investigate  the  transient  and  steady  state  flow  characteristics  of  an  HDR 
(Hot  Dry  Rock)  reservoir.   Each  of  the  problems  presented  here  use  the 
same  finite  element  mesh.  The  first  two  problems  are  steady  state  flow 
problems.   Fluid  is  not  allowed  to  leak  at  the  mesh  (reservoir)  boundary  for 
either  problem.   One  problem  has  the  injection  well  fluid  pressure  lower 
than  both  in-situ  reservoir  stresses;  the  other  problem  has  the  injection  well 
fluid  pressure  above  the  minimum  in-situ  reservoir  stress,  but  lower  than 
the  maximum  in-situ  reservoir  stress.    The  last  problem  has  transient 
boundary  conditions  that  model  experiment  #2070  conducted  by  Los  Alamos 
National  Laboratory  at  Fenton  Hill,  New  Mexico.   It  applies  a  pressure 
history  at  the  injection  well  (EE-3A)  and  monitors  the  pressure  rise  at  the 
shut-in  extraction  well  (EE-2).   There  is  no  far  field  fluid  leakage. 

5.1       Low  Pressure  Steady  State  Problem 

Figure  5.1  shows  the  mesh  used  for  this  steady  state  problem  with  the 
boundary  conditions  labelled.   Problems  with  only  pressure  boundary 
conditions  converge  more  rapidly  than  those  with  other  boundary  condition 
types.   For  this  problem,  with  two  pressure  boundary  conditions,  both  below 
the  minimum  in-situ  stress,  none  of  the  joints  will  completely  open  and  the 
joint  law  will  remain  in  effect. 
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Figure  5.2  is  a  fluid  pressure  contour  plot  showing  the  steady  state 
pressure  solution.   The  pressure  plot  is  almost  perfectly  symmetric  about  a 
vertical  center  line  through  the  mesh  with  little  pressure  drop  in  the  X 
direction  and  a  nearly  linear  pressure  drop  in  the  Y  direction.   The  flow 
rates  through  the  mesh  are  shown  in  Figure  5.3.   The  arrows  represent 
flow  direction  when  shown  on  the  plot.   The  arrows  are  not  shown  for  flow 
rates  out  of  the  legend  range.  The  flow  rates  are  constant  and  the  same  for 
each  of  the  three  vertical  flow  paths.   The  extremely  small  flow  rates  on  the 
boundaries  indicate  that  the  fluid  solution  is  converged  reasonably  well  for 
the  current  structure  solution.  A  careful  investigation  of  the  flow  rates  at 
the  extraction  well  shows  that  some  flow  does  move  around  behind  the 
extraction  well  to  exit  the  mesh  (see  Figure  5.4). 

Figure  5.5  is  a  displaced  mesh  plot  for  the  steady  state  flow  solution. 
Notice  that  the  higher  fluid  pressures  at  the  bottom  of  the  mesh  near  the 
injection  well  cause  both  horizontal  and  vertical  joints  to  open  more  than  at 
the  top  of  the  mesh.  Also  notice  that  the  stiffer  vertical  joints  are  not  as  far 
open  as  the  softer  horizontal  joints,  as  predicted  by  the  difference  in  the 
joint  opening  laws  for  the  two.   The  pressure  drop  is  low  enough  across  the 
soft  tensile  joints  that  they  provide  a  constant  pressure  "header"  across  the 
bottom  of  the  mesh. 

Throughout  the  mesh  the  displacements  are  not  as  uniform  as  the 
pressure  plot  might  indicate.   This  is  because  the  problem  is  not  completely 
converged.    To  reduce  computer  run  times  a  loose  structural  solution 
tolerance  was  used.   An  inspection  of  the  displacement  history  plots  for 
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several  nodes  in  the  mesh  shows  that  run  times  were  reduced  at  the 
expense  of  totally  converged  solutions.   Figure  5.6  compares  X  and  Y 
displacement  history  plots  for  three  nodes  in  the  mesh.   Each  of  the  X 
displacement  plots  show  good  convergence,  as  expected  from  the  displaced 
mesh  plot.   The  Y  displacement  history  plots  do  not  show  good  convergence 
even  though  the  specified  tolerance  was  met.   However,  the  solution  is 
nearly  converged.   The  line  plot  in  Figure  5.7  shows  the  Y  stress  varies  only 
slightly  in  the  mesh  from  the  top  to  the  bottom  and  is  near  the  expected 
value.   As  discussed  in  Section  3.1,  the  relatively  slow  convergence  is  due  to 
the  extremely  stiff  joint  law  that  forced  a  time  step  reduction  of  100.  A 
solution  to  this  difficulty  is  discussed  in  Chapter  VI. 

5.2       High  Pressure  Steady  State  Problem 

This  problem  has  specified  pressure  boundary  conditions  in  the  same 
locations  as  the  previous  steady  state  problem.   For  this  problem  the 
injection  pressure  is  15  MPa  and  the  extraction  pressure  is  still  1  MPa. 
When  the  fluid  pressure  rises  above  the  initial  in-situ  stress  the  joint 
elements  will  open  and  not  carry  any  load.   Only  the  fluid  pressure  in  the 
joints  will  act  on  the  rock  blocks,  which  can  "float"  in  the  middle  of  the 
mesh.   The  blocks  are  free  to  move  because  with  only  the  fluid  pressure 
acting,  translations  are  possible  without  changes  in  the  force  balance.   As 
fluid  pressures  continue  to  rise  and  the  joint  openings  increase,  the  blocks 
are  free  to  translate  more. 

In  this  problem  we  demonstrate  the  ability  to  model  floating  blocks  by 
specifying  in  inlet  pressure  of  15  MPa  and  an  outlet  pressure  of  1  MPa. 
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Because  the  inlet  pressure  is  above  the  minimum  in-situ  stress  of  10  MPa 
in  the  Y  direction,  joints  will  completely  open  and  blocks  will  "float"  in  the 
Y  direction  near  the  inlet. 

Figure  5.8  shows  the  calculated  pressure  distribution.    A  uniform 
pressure  does  exist  locally  around  the  inlet  location.   This  is  because  the 
joints  have  completely  opened  and  all  load  is  being  carried  by  the  fluid  in 
the  joints.   As  we  approach  the  outlet,  the  pressure  decreases  and  the 
structural  joints  are  carrying  part  of  the  load. 

To  satisfy  equilibrium,  the  final  stress  in  the  Y  direction  should  be  15 
MPa,  consistent  with  the  specified  input  pressure.   Figure  5.9  is  a  line  plot 
of  the  Y  stress  in  the  mesh.   The  Y  stress  is  not  constant  in  the  mesh  as 
expected.  Figure  5.10  is  a  line  plot  of  the  X  stress  and  shows  that  the  stress 
is  constant  along  the  line.   The  Y  stress  is  not  constant  along  the  line  of 
Figure  5.9  because  the  solution  is  not  completely  converged  due  to  a  loose 
convergence  tolerance.   As  the  solver  algorithm  seeks  a  solution,  the 
stresses  reach  local  equilibrium  with  the  fluid  pressures.    Global 
equilibrium  of  the  blocks  near  the  extraction  well  is  not  satisfied,  and  these 
blocks  must  translate  to  drive  the  boundaries  of  the  reservoir  out  against 
the  boundary  springs.   This  will  increase  the  loads  in  the  boundary  springs 
and  result  in  global  equilibrium.   Figure  5.9  shows  that  the  local  stresses  in 
the  rock  masses  have  come  up  to  the  current  fluid  solution,  but  the  rock 
masses  have  not  translated  to  bring  the  forces  and  stresses  to  equilibrium 
in  the  mesh.   The  Y  displacement  history  plots  in  Figure  5.11  show  this. 
Node  450  is  in  the  bottom  third  of  the  finite  element  mesh  and  node  1950  is 
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in  the  top  third  of  the  mesh.   Both  nodes  have  increasing  positive 
displacements,  but  from  Figure  3.4,  interface  material  three,  we  see  that 
the  displacements  at  the  top  of  the  mesh  need  to  be  more  than  twice  the 
current  displacement  of  node  1950  for  the  mesh  to  be  in  equilibrium. 

Again,  due  to  the  extremely  stiff  joint  model,  the  time  step  was 
reduced  by  a  factor  of  100.   As  a  result,  problem  run  times  for  complete 
convergence  would  be  excessive.   None  the  less,  it  appears  that  the  local 
solution  for  the  pressure  is  approximately  correct,  since  local  convergence 
has  been  obtained  and  the  translations  shown  in  Figure  5.10  are  in  the 
expected  direction. 

Two  factors  must  be  addressed  to  obtain  a  completely  converged 
solution.   First,  we  can  speed  the  solution  by  using  an  effective  joint  law 
with  the  same  flow  characteristics  but  a  softer  effective  structural  stiffness 
(Chapter  VI).   Secondly,  we  must  consider  whether  it  is  realistic  to  allow 
blocks  to  "float".   For  the  conditions  discussed  in  Section  5.2  the  blocks 
should  be  held  in  place  by  shear  forces  applied  by  the  X  normal  stress  that 
is  still  acting  on  the  blocks,  even  though  a  shear  model  is  not  implemented 
at  this  time.   It  may  be  more  realistic  to  tie  each  block  to  ground  to  limit 
rigid  body  motion.   Even  if  the  injection  pressure  is  above  both  in-situ 
stresses,  one  would  not  expect  blocks  as  large  as  the  ones  in  this  model  to  be 
completely  floating. 
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5.3       Transient  Flow  Problem 

This  problem  uses  the  same  finite  element  mesh  as  the  two  previous 
problems.   The  same  joint  opening  law  and  boundary  spring  stiffnesses  are 
also  used.   No  far  field  fluid  loss  is  allowed  at  the  reservoir  boundary.   As 
mentioned  before,  this  problem  models  a  constant  flow  rate  shut-in  test 
which  was  a  part  of  experiment  #2070  conducted  by  LANL  at  Fenton  Hill, 
New  Mexico.  Water  was  pumped  into  the  HDR  reservoir  at  a  flow  rate  of 
10.1e-3  m3/s  (133  gpm).  The  pressure  was  monitored  at  both  the  injection 
well  and  the  shut-in  extraction  well.   Figure  5.12  shows  the  measured 
pressure  history  at  the  injection  well  during  experiment  #2070. 

This  test  is  analogous  to  verification  problem  two,  where  a  single 
joint  was  pumped  up  like  a  balloon,  then  allowed  to  collapse.   For  the 
current  problem  the  shut-in  extraction  well  pressure  is  monitored  as  the 
reservoir  is  being  inflated.   Because  a  single  flow  rate  boundary  condition 
will  not  yield  a  unique  pressure  solution,  the  pressure  history  monitored  at 
the  injection  well  during  test  #2070  was  applied  and  the  resulting  input 
flow  rate  and  shut-in  extraction  pressure  were  monitored.   Figure  5.13  is  a 
plot  of  the  applied  injection  well  pressure  history  and  a  plot  showing  the 
application  point.   For  this  problem  the  bottom  of  the  mesh  is  considered  a 
line  of  symmetry  with  no  fluid  flow  through  this  boundary.   This  allows  us 
to  model  a  reservoir  that  is  twice  as  large  as  could  be  modelled  otherwise. 
Because  solution  times  are  basically  linear  with  problem  size,  the  run  time 
was  reduced  by  about  half. 
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Figure  5.14  is  a  pressure  contour  plot  for  the  first  time  step  of  the 
transient  fluid  solution.   Recall  that  the  injection  well  was  located  in  the 
lower  left  corner  of  the  mesh  and  that  no  flow  is  allowed  out  of  the  mesh. 
The  pressure  at  the  injection  well  did  rise  to  the  specified  pressure  and 
there  is  a  significant  pressure  gradient  from  that  point  to  the  top  of  the 
mesh.   The  minimum  pressure  in  the  mesh  is  0.97  MPa,  only  slightly  above 
the  initial  pressure  of  0.7  MPa. 

Figure  5.15  is  a  plot  comparing  the  measured  output  pressure  history 
for  well  EE-2  (shut  in  pressure  history)  to  the  calculated  pressure  history 
from  the  transient  boundary  conditions  discussed  above.   The  far  field  flow 
loss  was  zero,  so  all  the  flow  into  the  reservoir  was  accommodated  by  the 
joint  openings.   The  results  show  that  the  calculated  output  pressure 
increases  faster  than  the  measured  response. 

This  problem  is  an  initial  attempt  to  simulate  a  reservoir 
experiment.   It  demonstrates  that  the  features  needed  to  model  the 
experiment  are  working,  however,  a  careful  review  of  input  assumptions  is 
needed  to  develop  the  final  reservoir  model.   Three  factors  can  easily  slow 
the  transient  response  of  our  model:  increasing  the  reservoir  volume, 
introducing  a  far  field  flow  loss,  and  changing  the  time  step  size. 

Increasing  the  number  of  flow  paths  will  slow  the  transient  response 
by  producing  additional  fluid  storage  volume.  Experimental  data  indicates 
that  during  hydrofracturing  tests  on  the  Phase  II  reservoir,  99%  or  more  of 
the  injected  fluid  volume  was  accommodated  by  aseismic  tensile  fracturing 
(Brown,  1988b).  The  present  analysis  corroborates  this  observation,  with 
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nearly  all  input  flow  used  to  open  tensile  fractures.   If  more  tensile  joints 
are  included  in  the  model,  then  more  fluid  can  be  stored  in  the  tensile 
joints,  and  it  will  take  longer  for  the  fluid  to  reach  the  shut  in  well.   This 
will  slow  the  shut-in  transient. 

Introducing  a  far  field  loss  will  slow  the  transient  by  producing 
leakage  paths  for  the  fluid.  The  far  field  flow  loss  is  not  well  understood, 
but  researchers  at  LANL  (Brown,  1988a)  feel  that  the  far  field  flow  loss  for 
test  #2070  was  negligible  (less  than  2  %  of  the  injection  flow  rate).   In  the 
present  analysis,  no  far  field  loss  was  assumed. 

Finally,  changing  the  time  step  size  may  affect  the  response.   If  the 
time  step  in  the  solution  is  too  large,  the  solution  approaches  the  steady 
state  solution.  As  the  time  step  is  reduced,  the  results  converge  to  the 
transient  solution.   Additional  calculations  are  needed  in  which  we  vary 
the  time  step  size  to  examine  convergence. 

5.4       Summary 

Three  applications  problems  were  presented  in  this  chapter.   The  low 
pressure  steady  state  problem  showed  that  although  good  results  can  be 
obtained  for  a  large  problem,  care  must  be  taken  to  insure  that  the  problem 
is  converged.   The  high  pressure  steady  state  problem  confirmed  this.   With 
the  injection  well  pressure  above  the  in-situ  Y  stress,  the  blocks  were 
"floating"   in  the  Y  direction.   Large  translations  in  the  mesh  are  required 
of  most  of  the  mesh  to  compress  the  boundary  springs  enough  to  balance 
the  stress  applied  by  the  high  fluid  pressure.   However,  the  fluid  solution 
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should  not  be  significantly  affected  by  these  translations  because  the  mesh 
has  come  to  a  local  equilibrium  and  because  the  rigid  body  translation  does 
not  affect  the  three  major  shear  flow  paths  where  the  majority  of  the 
pressure  drop  occurs  between  the  inlet  an  outlet  wells.   The  transient 
problem  displayed  many  of  the  same  convergence  characteristics  as  the 
high  pressure  steady  state  problem  because  of  pressures  above  the  Y  in-situ 
stress.   The  transient  results  are  encouraging,  but  the  effects  of  the  fluid 
time  step,  reservoir  volume  and  far  field  flow  losses  need  to  be  investigated. 
The  solution  run  times  for  each  of  these  problems  was  significantly  slowed 
by  the  stiff  interface  elements. 
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P  =  1.0MPa 


P  =  7.5  MPa 


Figure  5.1:       Boundary  Conditions 

(low  pressure  steady  state  problem) 
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Figure  5.2:       Pressure  Contour  Plot 

(low  pressure  steady  state  problem) 
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Figure  5.3:       Flow  Rate  Contour  Plot 

(low  pressure  steady  state  problem) 


Figure  5.4:       Flow  Rate  Contour,  Zoom  Around  Extraction  Well 
(low  pressure  steady  state  problem) 
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Figure  5.5:       Displaced  Mesh  Plot 

(low  pressure  steady  state  problem) 


79 


NOOE    NUMBER    450 


-0.200-L 

0.000       0.500       1.000       1.500       2.000       2.500 


TIME 


(xlO  3) 


i 


NOOE  NUMBER  450 


0.100 


0.000 


1 


«  -0.200  + 


TIME 


Figure  5.6:       Displacement  History  Plots 

(low  pressure  steady  state  problem) 
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Figure  5.6  cont:    Displacement  History  Plots 

(low  pressure  steady  state  problem) 
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Figure  5.6  cont:    Displacement  History  Plots 

(low  pressure  steady  state  problem) 
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Figure  5.7a:     Y  Stress  Line  Plot 

(low  pressure  steady  state  problem) 
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Figure  5.7b:     Line  of  Application  for  Above  Line  Plot 
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Figure  5.8:       Pressure  Contour  Plot 

(high  pressure  steady  state  problem) 
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Figure  5.9a:     Y  Stress  Line  Plot 

(high  pressure  steady  state  problem) 
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Figure  5.9b:     Line  of  Application  for  Above  Line  Plot 


85 


TIME-   O.OOOOOE+OO 

-) 1 1- 


-3.000 -<• 

0.000  0.500  1.000  1.500  2.000  2.500  3.000 


(DISTANCE 


(xlO  2) 


Figure  5.10:     X  Stress  Line  Plot 

(high  pressure  steady  state  problem) 
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Figure  5.11:     Displacement  History  Plots 

(high  pressure  steady  state  problem) 
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Figure  5.13: 


b)  Application  Point 

Applied  Pressure  Boundary  Condition 
(transient  application  problem) 

89 


Figure  5.14: 


Pressure  Contour  Plot,  First  Transient  Step 
(transient  application  problem) 
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Figure  5.15: 


Comparison  of  Calculated  and  Experimental 
Output  Pressure  Histories 
(transient  application  problem) 
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Figure  5.15  ront:    Comparison  of  Calculated  and  Experimental 

Output  Pressure  Histories 
(transient  application  problem) 
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Chapter  VI 
Summary  and  Conclusions 

This  chapter  gives  a  summary  of  the  thesis  and  explains  the 
conclusions  that  can  be  drawn  from  the  verification  and  application 
problems.    Finally,  recommendations  are  made  for  future  work  both  for 
extending  DRACULA's  capabilities  and  for  modelling  HDR  reservoirs  with 
DRACULA. 

6.1       Summary 

In  Chapter  I  we  introduced  the  Hot  Dry  Rock  concept.   A  brief  history 
of  the  Hot  Dry  Rock  program  was  presented  to  better  understand  the  past 
modelling  and  experimental  results. 

In  Chapter  II  the  theoretical  basis  for  the  finite  element  method 
solution  was  developed  for  both  the  structural  and  fluid  models.   The 
structure  is  assumed  linear  elastic.   A  six  noded  quadratic  triangle  is  used 
to  approximate  the  structural  solution.   The  fluid  flow  rate  is  assumed 
proportional  to  the  pressure  gradient  and  the  joint  permeability,  which  is  a 
cubic  function  of  the  joint  openings.   A  three  noded  quadratic  line  element 
is  used  to  approximate  the  fluid  solution.   A  surface  element  is  presented 
which  when  superimposed  on  top  of  the  fluid  element  models  the  joint 
opening  law  influenced  by  fluid  pressures.   The  highly  nonlinear  nature  of 
the  problem  was  discussed  and  dynamic  relaxation  was  presented  as  a 
solution  method  capable  of  solving  the  problem. 
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Chapter  III  presents  the  computer  code  DRACULA  written  to 
implement  the  concepts  introduced  in  Chapter  II.   The  interactive  nature  of 
DRACULA  gives  the  user  the  ability  to  specify  a  problem  and  monitor  the 
results  all  in  the  same  interactive  graphics  environment.   Next,  some  input 
data  for  DRACULA  was  discussed  as  it  relates  to  the  boundary  conditions 
for  the  reservoir. 

In  Chapter  IV  we  presented  verification  problems  which  illustrate 
that  DRACULA  can  solve  highly  nonlinear  steady  state  and  transient  fluid 
flow  problems.  We  showed  that  flow  rate  and  far  field  flow  loss  boundary 
conditions  can  be  combined  in  one  problem  and  also  showed  that  the  choice 
of  the  time  step  size  for  transient  problems  is  important. 

Three  applications  of  DRACULA  were  illustrated  in  Chapter  V.   The 
first  problem  showed  that  good,  converged  results  on  these  nonlinear 
problems  are  attainable  even  on  large  problems.   Because  all  the  fluid 
pressures  were  below  the  in-situ  stress,  the  joint  laws  acted  as  expected 
with  softer  tensile  joints  opening  more  than  the  stiffer  shear  joints.   Very 
little  pressure  drop  was  seen  in  the  tensile  joints  compared  to  the  shear 
joints. 

The  second  and  third  problems  both  had  fluid  pressures  above  the 
minimum  in-situ  stress,  but  below  the  maximum  in-situ  stress. 
Convergence  was  reasonably  good  locally  where  fluid  pressures  were  above 
the  minimum  in-situ  stress,  but  rigid  body  translation  of  many  of  the  blocks 
in  the  mesh  are  required  before  global  convergence  can  be  attained. 
Excessive  run  times  are  required  for  these  problems  because  the  joint 
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stiffnesses  of  the  interface  elements  required  the  integration  time  step  to  be 
reduced  by  a  factor  of  100  to  maintain  stability.  The  third  problem  also 
showed  that  the  fluid  time  step  size,  reservoir  far  field  boundary  conditions 
and  number  of  fluid  flow  paths  should  be  evaluated  to  obtain  a  model 
transient  response  that  is  closer  to  the  experimental  transient  response. 

6.2       Conclusions 

An  analysis  tool  was  constructed  that  is  useful  for  analyzing  Hot  Dry 
Rock  reservoirs.   However,  a  careful  review  of  the  initial  and  boundary 
conditions  is  needed  to  develop  a  final  reservoir  model. 

Low  pressure  problems  are  easier  to  model  than  high  pressure 
problems  because  for  fluid  pressures  less  than  the  in-situ  stresses  the  joints 
do  not  open.   Therefore  the  mesh  is  not  required  to  translate  large  distances 
to  bring  the  mesh  to  equilibrium. 

High  pressure  problems  (those  with  fluid  pressures  above  the 
minimum  in-situ  stress)  converge  reasonably  well  at  a  local  level,  with 
local  structural  stresses  rising  to  local  fluid  pressures,  but  the  rigid  body 
modes  required  to  achieve  global  equilibrium  need  longer  run  times. 

The  transient  problem  was  essentially  a  set  of  high  pressure 
problems  and  displayed  the  characteristics  discussed  above.   For  the  mesh 
and  t-step  used,  the  transient  response  was  much  too  fast  compared  to  the 
experimental  data.   The  second  verification  problem  showed  that  reducing 
t-step  can  improve  a  solution.  The  effects  of  t-step  on  this  problem  should  be 
investigated. 
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The  tensile  joints  have  a  strong  effect  on  the  flow.  They  have  a  larger 
initial  opening  (slightly)  and  they  open  up  farther.   They  have  higher  flow 
rates  with  less  pressure  drop.   For  each  problem  investigated,  the  injection 
point  was  offset  in  the  mesh  to  attempt  to  force  flow  across  the  mesh,  but  for 
each  case  the  injection  point  location  seemed  to  have  little  bearing  on  the 
results.   The  pressure  horizontally  across  the  bottom  of  the  mesh  was 
essentially  constant  and  the  flow  rates  in  the  shear  paths  from  the  injection 
well  to  the  extraction  well  were  constant  and  equal  for  each  flow  path. 

Loose  convergence  tolerances  on  the  structure  gave  adequate  fluid 
results  for  the  low  pressure  problem,  but  the  fluid  results  for  the  high 
pressure  problem  may  change  slightly  when  the  problem  is  truly 
converged.   The  pressures  near  the  extraction  well  are  fixed  at  1.0  MPa. 
The  joint  stress  must  increase  to  about  14  MPa  for  a  total  stress  of  15.0  MPa 
to  balance  the  injection  fluid  pressure.   Stresses  this  high  in  an  interface 
element  will  require  joint  closure.  This  will  be  true  for  both  the  tensile  and 
shear  joints.   These  smaller  displacements  will  definitely  affect  the  flow 
solution,  probably  by  causing  higher  pressure  gradients  near  the  extraction 
well. 

The  joint  law  is  =100  times  stiffer  than  the  structural  elements.   This 
necessitates  reducing  the  time  step  by  a  factor  of  100  to  insure  numeric 
stability.   If  the  spring  stiffness  could  be  reduced  without  affecting  the  fluid 
solution  the  time  step  reduction  factor  could  be  increased  and  solution 
times  significantly  reduced. 
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Even  though  the  tolerance  for  a  problem  may  be  reasonable  for  a  good 
solution  for  some  problems,  other  problems  may  require  tighter  tolerances 
for  good  solutions.   The  results  from  the  two  steady  state  problems  show 
this.   The  problem  with  the  7.5  MPa  injection  pressure  converged  to  a 
reasonable  solution  with  a  structure  tolerance  of  O.le-3  and  a  fluid  tolerance 
of  0.5e-4.   However,  the  15  MPa  injection  pressure  problem  did  not  converge 
with  the  same  tolerances  and  the  same  damping  and  coupling  parameters. 
This  problem  required  tighter  tolerances  and  less  damping  for  better 
convergence. 

6.3       Recommendations 

As  discussed  in  Section  3.1  the  solution  algorithm  is  significantly 
slowed  by  the  extremely  stiff  joint  opening  law.   One  method  of  removing 
this  constraint  is  to  artificially  soften  the  joint  stiffness  for  the  global 
structural  solution,  but  still  recover  the  joint  openings  using  the  original 
joint  opening  law.   This  could  be  implemented  by  putting  a  relatively  soft 
spring  in  series  with  the  joint  law  for  the  structural  solution.   This  would 
not  significantly  affect  the  structural  solution.   If  this  method  were  to  work 
the  solution  time  could  be  reduced  by  a  factor  of  100. 

The  results  of  the  transient  problem  in  Section  5.3  are  encouraging 
from  the  stand  point  that  transient  results  were  obtained  for  such  a  large 
highly  nonlinear  problem,  but  the  results  did  not  match  well  with  the 
experimental  results  from  test  #2070  conducted  by  LANL.   Three  major 
factors  can  affect  the  transient  results  of  these  transient  problems:   (1 )  the 
fluid  time  step  size,  (2)  the  reservoir  volume  and  the  distribution  of  volume 
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between  shear  and  tensile  fractures,  and  (3)  far  field  fluid  boundary 
conditions.   Each  of  these  should  be  studied  to  better  understand  their 
effects  on  transient  problems. 

The  results  of  this  thesis  are  a  foundation  on  which  to  build  a  fully 
coupled  heat  transfer,  fluid  flow,  rock  deformation  analysis  code.   The  fluid 
flow  and  rock  deformation  models  are  now  fully  coupled.   The  next  step  is  to 
include  temperature  dependent  fluid  viscosity,  thermal  stresses  and  a  heat 
transfer  model  for  heat  flow  between  the  structure  and  the  fluid. 

Dynamic  relaxation  has  proven  to  be  a  very  reliable  solution  method 
for  these  nonlinear  problems.   In  addition  it  has  the  virtue  of  approximately 
linear  solution  time  increase  with  problem  size.   However,  it  is  not  a  fast 
solution  method.   Many  other  methods  exist  that  are  reasonably  robust  with 
good  initial  guesses  on  the  solution  and  are  much  faster  than  dynamic 
relaxation.    Since  dynamic  relaxation  seems  to  be  at  its  worst  when  it  is 
near  a  solution,  other  methods  should  be  investigated  with  the  intent  that 
dynamic  relaxation  would  start  the  problem  and  an  alternate  solution  used 
for  final  convergence.  The  effects  of  the  rigid  body  modes  of  many  of  the 
blocks  in  the  finite  element  mesh  must  be  considered  when  investigating  a 
new  solution  algorithm. 

Some  of  the  "floating"  block  problems  described  in  Section  5.3  would 
not  present  themselves  if  a  shear  law  existed  in  the  model.   If  a  shear 
model  is  implemented,  DRACULA  should  be  modified  such  that  the  joint 
laws  are  entered  in  terms  of  the  principal  stress  orientation  rather  than  as 
data  in  a  table.   This  would  allow  the  effects  of  the  principal  stress 
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orientation  to  be  investigated  much  more  easily  than  with  the  current  joint 
law  implementation. 
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ABSTRACT 

Hot  Dry  Rock  geothermal  reservoirs  are  currently  being  developed 
around  the  world.   Extensive  experimental  work  has  been  done,  but  few 
numerical  models  exist.   This  model  solves  the  coupled  fluid  flow  --  rock 
deformation  problem  and  provides  the  foundation  for  fully  coupled  heat 
transfer  with  thermal  stresses. 

A  computer  program  called  DRACULA  (Dry  Rock  Analysis  Code)  is 
used  to  model  steady  state  and  transient  fluid  flow  results.   The  model  uses 
discrete  joints  to  model  the  fluid  flow  through  the  fractured  rock.   Two 
separate  models  are  used  in  the  program:  (1 )  a  structural  model  to  solve  for 
rock  deformations,  and  (2)  a  fluid  model  to  solve  for  fluid  pressures.   The 
two  models  are  coupled  through  the  joint  permeability  (cubic  law)  and  the 
effective  stress  law.   Dynamic  relaxation  is  used  to  obtain  solutions. 

The  formulation  includes  the  transient  storage  terms  associated  with 
joint  opening  velocities.   This  allows  the  user  to  simulate  reservoir 
operation  in  an  inflation/deflation  mode  where  fluid  is  pumped  into  the 
reservoir,  stored  by  opening  rock  joints,  and  then  recovered  when  the 
pressure  is  lowered. 

Verification  problems  demonstrate  that  the  fluid  and  structural 
models  are  correctly  coupled.    The  first  verification  problem  shows  flow 
through  a  single  flow  path  where  joint  opening  is  not  uniform.   The  second 
illustrates  inflation/deflation  in  a  single  flow  path. 


Results  are  presented  for  a  more  complex  model  simulation  of  the 
Fenton  Hill  reservoir.   Results  for  flow  between  two  wells  at  specified 
pressures  show  that  the  tensile  joints  (normal  to  the  smallest  in-situ  stress) 
open  and  have  higher  flow  rates  with  less  pressure  drop  than  the  shear 
joints  (normal  to  the  maximum  in-situ  stress).   When  the  injection 
pressure  is  raised  above  the  minimum  in-situ  stress,  the  tensile  joints 
completely  open  and  rock  masses  "float".    An  initial  simulation  of  a  shut-in 
experiment  at  Fenton  Hill  demonstrates  the  storage  effect  of  joint  openings. 

Recommendations  are  made  on  speeding  solutions  and  further 
calculations  needed  to  model  the  Fenton  Hill  reservoir. 


